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ABSTRACT 


Viscous  effects  in  low-density  nozzle  flows  were  investigated 
numerically,  and  comparisons  were  made  with  experimental  data.  The 
numerical  method  of  Patankar  and  Spalding  was  modified  to  solve  the 
internal  laminar  bound  ary -layer  equations  for  two-dimensional  flow  or 
axisymmetric  flow  with  or  without  transverse  curvature.  A  listing  is 
given  of  the  computer  code.  Bound  ary -layer  displacement  thicknesses 
for  typical  nozzle  geometries  and  flow  conditions  are  presented. 
Solutions  were  obtained  for  specific  conditions  corresponding  to  experi¬ 
mental  data.  The  result  is  a  relatively  fast,  simple  to  use  numerical 
procedure,  which  is  shown  to  give  results  in  good  agreement  with 
experimental  data. 


111 


AEDC-TR-72-52 


CONTENTS 

Page 

ABSTRACT . iii 

NOMENCLATURE . vii 

I.  INTRODUCTION .  1 

II.  BASIC  EQUATIONS 

2.  1  Governing  Equations .  3 

2.  2  Nondimensionalized  Equations  . .  5 

2.3  Transformed  Equations .  7 

III.  NUMERICAL  SOLUTION 

3.  1  Basic  Scheme .  9 

3.2  Input  Conditions .  11 

3.3  Some  Possible  Computational  Difficulties .  14 

IV.  NUMERICAL  RESULTS .  16 

V.  COMPARISONS  WITH  EXPERIMENTAL  DATA  AND  OTHER 

THEORETICAL  RESULTS .  18 

VI.  CONCLUSIONS . 21 

REFERENCES . 22 

APPENDIXES 

I.  ILLUSTRATIONS 
Figure 

1.  Definition  of  Coordinate  System . 27 

2.  Displacement  Thickness  of  Successive  Iterations  for 

Condition  No.  5  (Table  I) . 28 

3.  Displacement  Thickness  for  Condition  No.  1 

(Table  I) . 29 

4.  Displacement  Thickness  for  Condition  No.  2 

(Table  I) . 30 

5.  Displacement  Thickness  for  Condition  No.  3 

(Table  I) . 31 

6.  Displacement  Thickness  for  Condition  No.  4 

(Table  I) . 32 


v 


AEDC-TR-73-52 


F  igure  Page 


7.  Displacement  Thickness  for  Condition  No.  5 

(Table  I) .  33 

8.  Displacement  Thickness  for  Condition  No.  6 

(Table  I) .  34 

9.  Displacement  Thickness  for  Condition  No.  7 

(Table  I) .  35 

10.  Displacement  Thickness  for  Condition  No.  8 

(Table  I) .  36 

11.  Nozzle  Exit  Velocity  and  Temperature  Profiles  for 

Condition  Nos.  3  and  9  (Table  I) .  37 

12.  Displacement  Thicknesses  of  an  Axisymmetric  and 

Two-Dimensional  Nozzle  with  Equal  Area  Ratios  ....  38 

13.  Reservoir  Reynolds  Number  as  a  Function  of  pQ,  TQ, 

and  r5>'  for  Various  Gases .  39 

14.  Calculated  and  Measured  Pitot  Pressure  Profiles  at  the 

M3  Nozzle  Throat  for  ReQj  r*  =  900  .  40 

15.  Calculated  and  Measured  Pitot  Pressure  Profiles  at  the 

M3  Nozzle  Exit  for  Re0j  r*  =  1800  .  41 

16.  Calculated  and  Measured  Pitot  Pressure  Profiles  at  the 

Nozzle  Exit  for  ReQ  r*  =  900  .  42 

17.  Calculated  and  Measured  Pitot  Pressure  Profiles  at  the 

M3  Nozzle  Exit  for  ReQ  r*  =  4500  .  43 

18.  Calciilated  and  Measured  Pitot  Pressure  Profiles  Down¬ 
stream  of  the  M3  Nozzle  Throat  for  ReQ  r*  =  1170  ...  44 

19.  Calculated  and  Measured  Pitot  Pressure  Profiles  Down¬ 
stream  of  the  M3  Nozzle  Throat  for  Re0j  r*  =  330  ....  45 

20.  Calculated  Velocity,  Temperature,  and  Pitot  Pressure 
Profiles  Downstream  of  the  M3  Nozzle  Throat  for 

Re0<r*  =  330  . 46 

21.  Calculated  and  Measured  Heat-Transfer  Coefficients  in  a 

Low -Density  Nozzle .  47 


VI 


AEDC-TR-72-52 


Page 

II.  TABLE 

I.  Conditions  for  the  Numerical  Solutions  Presented  in 


Section  IV  .....  . .  48 

III.  BOUNDARY -LAYER  COMPUTER  CODE .  49 


IV.  INITIAL  PRESSURE  DISTRIBUTION  COMPUTER  CODE  .  .  67 

NOMENCLATURE 


A  , 

CP 

Cv 

d* 

H  : 

H 

h 

k  .. 

M 

• 

m 

Pr 


R 


Re  * 
neo,  r 


Cross-sectional  area  of  nozzle 
Specific  heat  at  constant  pressure 
Specific  heat  at  constant  volume 
Nozzle  throat  diameter 
Local  total  enthalpy 

h/hg 

Heat-transfer  coefficient 
Thermal  conductivity 
Mach  number 
Mass  flux 

Prandtl  number,  pcp/k 
Pressure 

Measured  pitot  pressure 
(7  -  Dp/  (7p0) 

Nozzle  wall  heat -transfer  rate 

Longitudinal  radius  of  curvature  of  converging  portion  of 
nozzle 

R/rT 

Nozzle  reservoir  Reynolds  number,  p0(2Ho)  r*/M0 


vii 


AEDC-TR-72-52 


r 

r 

r"' 


T 

T 


aw 


T 


u 

u 

V 

V 
X 
X 


y 

y 

z 

z 

or 

7 

6 


6" 

c 

>1 

P 

P 

1/ 

I 

P 

P 


Defined  by  Eq.  (4) 

/  * 
r/ r 

Nozzle  throat  radius 
Temperature 

Adiabatic  wall  temperature 
T/T0 

Velocity  component  in  x  direction 
u/(Hq)1/2 

Velocity  component  in  y  direction 
v/(Hq)1/2 

Coordinate  along  nozzle  wall 

/  * 
x/r 

Coordinate  normal  to  nozzle  wall 

y/r 

Coordinate  along  nozzle  axis  referenced  to  the  throat 
z/r"' 

Nozzle  wall  angle 

Cp/Cy 

Boundary  -layer  thickness  defined  as  the  value  of  y  where 
u/uE  =  0.  99 

Boundary-layer  displacement  thickness 
Exponent  in  power  law  viscosity,  p  ~ 

y/yE 

Viscosity 

P/P0 

Denotes  planar  or  axisymmetric  flow  in  Eq.  (23) 
Transformed  x  coordinate,  Eq.  (11) 

Mass  density 

P/P0 


viii 


AEDC-TR-73-52 


ip  Stream  function 

u  Transformed  stream  function,  Eq.  (12) 


SUPERSCRIPT 

Condition  immediately  downstream  of  a  normal  shock 


SUBSCRIPTS 


<L 

E 

e 

I 

o 

w 


Nozzle  centerline 

Outer  edge  of  boundary  layer 

Nozzle  exit 

Inner  edge  of  boundary  layer 
Reservoir  (total)  conditions 
Nozzle  wall 


IX 


AEDC-TR-73-52 


SECTION  I 
INTRODUCTION 


Investigations  of  viscous  effects  in  low -density  gas  flows  in  two- 
dimensional  and  axisymmetric  channels  and  nozzles  have  been  con¬ 
ducted  during  the  past  few  years  in  support  of  the  design  of  low -density 
wind  tunnel  facilities  and  small  microthrust  rockets  used  for  space¬ 
craft  attitude  control.  In  addition  to  these  areas  of  interest,  there  are 
two  other  areas  where  attention  to  viscous  effects  is  required.  One  is 
that  of  internal  bound  ary -layer  scaling.  When  the  nozzle  and/or  plume 
flow  of  a  rocket  engine  is  to  be  investigated  in  a  wind  tunnel  or  space 
chamber,  it  is  usually  necessary  to  significantly  reduce  the  size  of  the 
nozzle  used.  To  achieve  adequate  simulation,  the  model  nozzle  viscous 
effects  must  appropriately  simulate  those  of  the  actual  rocket.  A  spe¬ 
cific  example  of  this  problem  is  found  in  the  laboratory  study  of  rocket 
exhaust  plumes  interacting  with  the  free -stream  (Ref.  1).  The  other 
area  of  interest  is  associated  with  the  study  of  gas  dynamic  and  chemi¬ 
cal  lasers  (Ref.  2).  The  total  laser  power  output  is  influenced  by  the 
static  gas  pressure  in  the  optical  cavity  just  downstream  of  the  nozzle 
exit.  Since  the  nozzles  used  are  small  and  since  they  operate  at  rela¬ 
tively  high  total  temperatures,  the  nozzle  viscous  effects  significantly 
influence  the  nozzle  static  pressure  distribution.  Although  numerous 
comparisons  with  experimental  data  from  low -density  wind  tunnel 
nozzles  will  be  presented  to  check  the  accuracy  of  the  results,  the 
present  investigation  was  motivated  by  the  current  interest  in  rocket 
nozzle  scaling  and  the  viscous  effects  on  the  operation  of  gas  dynamic 
and  chemical  laser  systems. 

Some  previous  investigations  of  the  design  and  analysis  of  low- 
density  wind  tunnel  nozzles  are  given  in  Refs.  3  through  5.  The 
method  of  Potter  and  Carden  (Ref.  3)  was  developed  to  design  nozzles 
for  particular  test  section  flow  conditions.  It  is  based  on  an  integral 
technique  which  uses  the  similar  solutions  of  Cohen  and  Reshotko 
(Ref.  6).  Although  the  work  of  Potter  and  Carden  (Ref.  3)  has  proved 
successful  for  the  design  of  nozzles  to  produce  desired  flow  conditions, 
it  is  not  directly  applicable  to  the  analysis  of  specified  nozizle  geometry. 
Also,  non-similar  solutions  presented  in  Ref.  4  indicate  that  similarity 
does  not  exist  in  nozzle  flows,  particularly  near  the  throat,  and  as  a 
consequence,  inaccuracies  may  result  from  the  use  of  similar  solutions 
for  relatively  short  nozzles. 
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The  method  presented  in  Refs.  4  and  5  solves  the  non-similar 
laminar  boundary-layer  equations  with  or  without  first-order  transverse 
curvature  (referred  to  as  second-order  in  Refs.  4  and  7),  with  or  with¬ 
out  velocity  slip  and  temperature  jump  boundary  conditions,  and  it  has 
been  shown  to  be  accurate.  However,  the  method  has  certain  disadvan¬ 
tages:  (1)  the  numerical  integration  scheme  is  that  of  Jaffe,  Lind,  and 
Smith  (Ref.  7)  which  requires  a  relatively  large  amount  of  computer 
time,  (2)  the  computer  program  is  large  and  not  necessarily  simple  to 
use,  (3)  the  transformation  variables  are  not  amenable  to  internal  flow 
problems,  and  (4)  large  flow  expansions  frequently  require  interpo¬ 
lating  the  solutions,  changing  the  numerical  step  size  across  the  nozzle, 
and  resuming  the  calculations.  These  disadvantages  discourage  the  use 
of  the  method  described  in  Refs.  4  and  5. 

The  methods  described  in  Refs.  3  through  5  assume  that  the  nozzle 
flow  consists  of  a  viscous  region  and  an  inviscid  (core)  region.  Some 
previous  investigations  which  are  not  restricted  to  such  flows,  but  per¬ 
mit  viscous  effects  across  the  entire  channel,  are  described  in  Refs.  8 
through  10.  These  investigations  are  more  suitable  for  the  study  of 
flow  in  microthrust  rockets  where  the  flow  may  be  fully  viscous.  How¬ 
ever,  these  methods  also  have  certain  disadvantages  for  the  present 
application.  The  works  of  Adams  (Ref.  8)  and  Williams  (Ref.  9)  are 
based  on  similar  solutions,  the  conditions  for  which  are  not  likely  to  be 
satisfied  by  a  large  class  of  practical  problems.  The  work  of  Rae 
(Ref.  10)  is  a  significant  contribution  to  the  study  of  low -density  nozzle 
flows,  and  this  method  will  be  discussed  further  in  this  report. 
Numerical  results  from  Rae's  method  will  be  compared  with  results 
from  Ref.  5,  results  obtained  in  the  present  investigation,  and  experi¬ 
mental  data.  Rae’s  method  was  shown  to  give  good  results  for  fully 
viscous  flows  (Ref,  10);  however,  it  will  be  shown  herein  to  be  less 
accurate  than  the  method  of  Refs.  4  and  5  or  the  present  method  for  the 
flow  regime  of  interest  in  this  investigation. 

The  objectives  of  the  present  investigation  are:  (1)  to  provide 
results  for  estimating  the  viscous  effects  in  low-density  converging- 
diverging  nozzle  flows  based  on  certain  flow  parameters  and  nozzle 
geometries,  and  (2)  to  provide  a  fast,  simple  to  use  method  for  calcu¬ 
lating  viscous  effects  in  low-density  nozzle  flows.  To  meet  these  ob¬ 
jectives,  the  numerical  integration  scheme  of  Patankar  and  Spalding  . 
(Ref.  11)  was  used.  In  fact,  the  program  used  for  this  work  was  taken 
directly  from  Ref.  11  and  then  suitably  modified  for  internal,  low- 
density  flows.  The  resulting  program  is  similar  to  that  of  Mayne 
(Ref.  2)  except  for  the  numerical  scheme  used  near  the  wall.  Also,  the 
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present  work  is  based  on  a  set  of  dimensionless  equations  which  reduces 
the  amount  of  input  and  provides  more  convenient  solutions  in  the  sense 
that  they  are  applicable  to  flows  which  satisfy  certain  parameters 
rather  than  specific  inputs  of  pressure,  temperature,  etc. 

The  following  section  describes  the  governing  equations  and  boundary 
conditions,  the  variables  used  to  nondimensionlize  the  equations,  and  the 
transformed  equations.  Section  III  describes  briefly  the  numerical  solu¬ 
tion  as  well  as  some  of  the  computational  difficulties  which  have  been 
encountered  with  this  program.  Section  IV  is  devoted  to  numerical  re¬ 
sults  and  Section  V  to  comparisons  of  the  present  results  with  previous 
theoretical  investigations  and  experimental  data.  Some  conclusions  are 
given  in  Section  VI. 


SECTION  II 
BASIC  EQUATIONS 


The  governing  system  of  equations  and  the  boundary  conditions  are 
presented  in  this  section.  Certain  dimensionless  and  transformation 
variables  are  introduced  and  used  to  transform  the  governing  equations 
for  convenience  of  the  numerical  solution.  Also,  the  boundary-layer 
displacement  thickness  is  derived  in  the  transformed  plane.  Although 
the  program  will  solve  the  two-dimensional  equations  or  axisymmetric 
equations  with  or  without  transverse  curvature,  only  the  equations 
appropriate  to  the  latter  with  the  transverse  curvature  terms  retained 
are  considered  here  in  detail  because  they  are  more  general.  The  two- 
dimensional  equations  can  be  obtained  from  the  equations  considered  by 
setting  r  =  1,  and  the  axisymmetric  equations  without  transverse  curva¬ 
ture  can  be  obtained  by  setting  r  =  rw. 


2.1  GOVERNING  EQUATIONS 

The  governing  system  of  equations  is  taken  as  that  obtained  by 
Probstein  and  Elliott  (Ref.  12).  The  equations  in  curvilinear  coordinates 
are: 

Continuity  Equation 


d(pur)  d(pvr) 
<9x  +  dy 


(1) 
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Momentum  Equation 


(9 H  du  dp  1  <9  /  <9u\ 
+  pvTy  =  “77  +  7  57^57; 


Total  Energy  Equation 


dll  dll 

puT*  +  pvdV 


l  d_\  [fi  dll 
t  dy  Y  Pr  <9y 


(2) 


(3) 


The  total  energy  equation  is  obtained  by  multiplying  the  momentum 
equation  by  u  and  adding  the  result  to  the  static  energy  equation.  The 
coordinate  system  is  defined  in  Fig.  1  (Appendix  I)  with  the  r(x,  y) 
term  defined  for  internal  flow  as 

r(x,>)  -  rtt(x)  -  >  cos  a  ^ 


Probstein  and  Elliott  (Ref.  12)  obtained  Eqs.  (1)  through  (3)  by  an 
order  of  magnitude  analysis  of  the  general  forms  of  the  continuity  equa¬ 
tion,  Navier- Stokes  momentum  equations,  and  energy  equation.  The 
assumptions  made  in  the  analysis  were  that  the  ratio  of  the  boundary- 
layer  thickness  to  the  longitudinal  radius  of  curvature  of  the  body  sur¬ 
face  was  small  compared  to  unity,  and  the  ratio  of  the  boundary- 
layer  thickness  to  the  nozzle  radius  was  on  the  order  of  unity. 
Therefore,  Eqs.  (1)  through  (3)  are  valid  for  nozzles  which  have  a 
longitudinal  radius  of  curvature  much  larger  than  the  nozzle  radius. 

This  stipulation  is  normally  satisfied  by  axisymmetric  convergent - 
divergent  nozzles  used  for  low -density  wind  tunnels,  rockets,  and 
gas  dynamic  and  chemical  laser  systems. 

The  axisymmetric  boundary-layer  equations  without  transverse 
curvature  terms  correspond  to  those  which  can  be  obtained  from 
Eqs.  (1)  through  (3)  by  replacing  r(x,  y)  with  rw(x)  as  stated  above. 
Because  rw(x)  is  a  function  of  x  only,  it  can  be  eliminated  from  Eqs.  (2) 
and  (3),  and  therefore  appears  only  in  the  continuity  equation.  The  re¬ 
sulting  set  of  equations  can  be  used  to  describe  internal  or  external 
boundary  layers.  It  was  shown  in  Ref.  5,  by  solutions  with  and  without 
the  transverse  curvature  terms,  that  the  effect  of  transverse  curvature 
is  important  for  6/rw(x)  ~  0(1). 
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Implicit  in  Eq.  (2)  by  the  use  of  the  total  derivative  of  p  with  re¬ 
spect  to  x  is  the  y  component  of  the  momentum  equation  as  given  by 


The  validity  of  this  equation  is  sometimes  questioned  for  thick  laminar 
boundary  layers.  However,  the  analysis  of  Probstein  and  Elliott 
(Ref.  12)  indicates  that  this  equation  is  consistent  with  the  other  equa¬ 
tions  in  the  governing  set.  Equation  (5)  was  used  in  Refs.  4  and  5, 
and  good  agreement  between  calculated  and  measured  boundary -layer 
profiles  was  obtained  for  flow  conditions  where  99  percent  of  the  cross- 
sectional  area  of  a  nozzle  was  boundary  layer. 

For  the  boundary  conditions  at  the  edge  of  the  boundary  layer,  it 
was  assumed  that  an  isentropic  core  flow  exists  along  the  nozzle  center- 
line,  from  which  the  flow  properties  at  the  edge  of  the  boundary  layer 
can  be  calculated.  The  boundary  conditions  at  the  nozzle  wall  were' 
taken  as  zero  velocity  and  a  prescribed  wall  temperature  distribution 
or  a  prescribed  wall  heat-transfer  distribution.  Solutions  were  pre¬ 
sented  in  Ref.  5  with  and  without  velocity  slip  and  temperature  jump 
boundary  conditions.  For  the  conditions  investigated  in  Ref.  5,  it  was 
found  that  the  nozzle  flow  merged  (that  is,  the  boundary  layer  com¬ 
pletely  filled  the  nozzle)  before  velocity  slip  and  temperature  jump 
became  significant.  Therefore,  no-slip  wall  boundary  conditions  were 
assumed  adequate  for  this  investigation. 

The  system  of  equations  is  completed  by  using  the  equation  of  state, 
p  =  pRT,  and  expressing  the  viscosity  as  some  function  of  T,  p  =  p(T). 
The  governing  equations  are  next  made  dimensionless. 


2.2  NONDIMENSIONAL1ZED  EQUATIONS 

Dimensionless  variables  are  used  to  nondimens ionalize  the  govern 


ing  equations  as  follows: 


(6d) 

P  <&• 

(6g) 

(6e) 

“ '  W 

(6h) 

(6f ) 

V  "  lu  Vs 

(6i) 
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Using  these  equations  in  Eqs.  (1)  through  (3)  gives  for  the  continuity, 
momentum,  and  total  energy  equation,  respectively. 


and 


where 


(7) 

(8) 


(9) 


(10) 


A  few  comments  are  in  order  concerning  the  choice  of  the  dimen¬ 
sionless  variables.  Note  that  (Hq)1^2  was  used  to  normalize  the  velocity 
components  u  and  v  instead  of  the  more  commonly  used  maximum 
velocity  (2H0)*'  2.  The  motivation  for  this  was  to  recover  the  same 
form  of  the  boundary -layer  equations  in  dimensionless  variables  as  in 
physical  plane  variables.  If  (2H0)^2  were  used  in  Eqs.  (6h)  and  (6i)  in 
place  of  (H0)l/2,  then  the  term  M  [1  -  (1/Pr)]u  9u/ 8y  in  the  right  hand 
side  of  Eq.  (9)  would  be  2ju  [1  -  (l/Pr)]u9u/9y ,  which  would  not  be  the 
same  as  in  Eq.  (3).  The  maximum  velocity  (2H0)^/2  was  used  in  the 


This  introduces  the  term  (2) 


1/2 


definition  of  Re0  r*  for  convenience, 
as  a  coefficient  of  Re0)  r*,  but  since  Re0j  r*  is  a  constant  for  each 
solution,  the  constant  <2) “1/2  causes  no  inconveniences  in  the  numerical 
solutions.  Note  that  by  defining  a  new  viscosity,  P  =  (2)^'2P7Re0j  r*, 
Eqs.  (7)  through  (9)  are  identical  to  Eqs.  (1)  through  (3).  This,  in  fact, 
was  done  in  the  numerical  computations. 


For  a  given  set  of  initial  conditions,  boundary  conditions,  y,  Pr, 
nozzle  geometry,  etc. ,  the  solution  to  Eqs.  (7)  through  (9)  depends 
only  on  ReQ  r*  if  a  power  law  viscosity  variation  with  temperature  is 
assumed,  i/e.,  p  =  T*  =(T/T0)S.  However,  if,  for  example, 
Sutherland's  viscosity  law  is  assumed,  then  the  solutions  will  also 
depend  on  the  absolute  value  of  temperature.  Because  the  calculation  of 
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viscosity  is  carried  out  in  a  separate  subroutine  in  the  computer  pro¬ 
gram,  it  is  convenient  to  use  any  viscosity  law  desired.  Solutions  will 
be  presented  for  both  Sutherland's  law  and  power  law  variations  of 
viscosity. 


2.3  TRANSFORMED  EQUATIONS 


Equations  (7)  through  (9)  were  transformed  from  the  dimensionless 
physical  x  -  y  plane  to  the  §  -  u  plane  by  the  transformation  variables 

=  X  (11) 


and 


<d(x,y)  = 


ip(x, y)  -  li^tx) 

^eg<)  - 


(12) 


where  ^(x,  y)  is  the  stream  function  which  identically  satisfies  the  con¬ 
tinuity  equation,  Eq,  (7),  i.  e. , 


dip  _ 

T*  m  -pvr 


(13) 


and 


dtp  - 

57  -  pur 


(14) 


The  transformation  used  is  that  due  to  von  Mises  with  u  introduced  to 
restrict  the  integration  across  the  boundary  layer  from  zero  to  unity. 
The  subscripts  I  and  E  denote  the  inner  and  outer  edge  of  the  boundary 
layer,  respectively,  and  fi  and  ^  are  functions  of  x  only.  The 
operators  (8/ 83c )  and  (8/ 8y)  are  given  by 


d_ 

87 


d 


do  d 
8x  do 


d 

3f 


dca 

Wi 


dip  l 

IT 


do 

dtp  v 


W)  d_ 

<fx  /  dot 


(15) 


and 


d  _  do  d  do  dip  d 

d~y  dy  do  dip  dy  do 


(16) 
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Using  Eqs.  (12)  through  (14)  in  Eqs.  (15)  and  (16)  and  applying  them  to 
Eqs.  (8)  and  (9)  give  for  the  momentum  and  total  energy  equations 


and 


dH  Vjh  +  KrE™E 


(2)' 


Re 


d_ 

da> 


_2 _ 

r  P  u  M 


1  dp 
pH  d"x 


dH 

rjiiij  +  to  (r^ny  —  r^m^) 

(911  (2)Ji 

~  _2 _  _T 

i  p  up  <9h 

cK  + 

dto  11  e  Q 

dto 

Pr  doi 

where 


+ 


d 

2 

T  p  U  p 

dto 

J<£E  -  W,)2 

d"x 


-  rjirij 


and 


d^E 

IT”  =  “rEmK 


(17) 


(18) 


(19) 


(20) 


and  where  £(x)  has  been  replaced  by  x.  Except  for  the  constant 

(2)^2/Re0> r*,  Eqs.  (17)  and  (18)  are  identical  to  the  momentum  and 
total  energy  equations  for  external  flow  which  are  solved  numerically 
by  the  method  of  Patarikar  and  Spalding  (Ref.  11).  The  calculation  of 
the  mass  transfer  fluxes  mj  and  mg  will  be  discussed  in  Section  III. 


The  boundary  -layer  displacement  thickness  was  used  to  calculate 
an  effective  nozzle  geometry  which  in  turn  was  used  to  calculate  an 
axial  pressure  distribution.  The  displacement  thickness,  6*,  which 
takes  into  account  transverse  curvature  is  expressed  by  Probstein  and 
Elliott  (Ref.  12)  as 


/  2OTpEuEdy  =  J"  2?rr(pEUp  -  pu)dy 

o  o 


(21) 
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Solving  Eq.  (21)  for  6  gives 


A  quadratic  equation  for  6  must  be  solved  to  obtain  Eq.  (22)  from 
Eq.  (21).  The  positive  sign  was  chosen  so  that  S*/r*  <  rw/cos  a. 
The  displacement  thicknesses  for  axisymmetric  flows  without  trans¬ 
verse  curvature  and  planar  flows  can  be  obtained  by  using  (rw)v  in 
place  of  r  in  Eq.  (21).  The  displacement  thicknesses  are  then  given 

by 


8*  _  Wz-tJ 

r*  ~  .v-  — 

(rw>  p EUE  (23) 

where  v  =  0  for  planar  flows  and  v  =  1  for  axisymmetric  flows  without 
transverse  curvature. 


SECTION  III 
NUMERICAL  SOLUTION 


The  numerical  solution  of  Eqs,  (17)  and  (18)  is  discussed  in  this 
section.  The  basic  scheme  is  briefly  described  and  a  listing  of  the 
computer  code  is  presented.  A  description  of  the  necessary  inputs  to 
the  program  is  given,  and  some  computational  difficulties  which  have 
been  encountered  with  the  program  are  discussed. 


3.1  BASIC  SCHEME 

As  previously  stated,  the  basic  numerical  integration  scheme  of 
Patankar  and  Spalding  (Ref.  11)  was  suitably  modified  and  used  to  solve 
the  governing  internal  flow  equations  pertinent  to  the  present  investiga¬ 
tion.  Because  the  flow  regime  of  the  present  investigation  was  entirely 
laminar,  some  subroutines  of  the  computer  code  listed  in  Ref.  11  which 
were  associated  with  turbulent  flow  were  removed  for  the  present  code. 
Also,  most,  but  not  all,  statements  pertaining  to  turbulent  flow  were 
removed.  The  computer  code  as  used  for  the  present  investigation  is 
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listed  in  Appendix  III.  Essentially  this  same  program  was  used  on  two 
different  computers  at  AEDC,  a  Scientific  Data  Systems  (SDS)  9300  and 
an  International  Business  Machines  (IBM)  370/155.  The  code  in 
Appendix  III  is  the  one  used  on  the  IBM  370/155. 

The  basic  numerical  scheme  used  by  Patankar  and  Spalding  is  dis¬ 
cussed  in  detail  in  Ref.  11.  One  of  the  primary  features  of  this  finite 
difference  technique  is  that  the  set  of  linear  algebraic  equations  which 
must  be  solved  has  only  three  unknowns  in  each  equation,  and  this  set 
can  be  solved  by  simple  successive  substitution  {Gauss  reduction,  or 
elimination.  Ref.  13)  rather  than  by  matrix  inversion.  This  technique 
provides  a  saving  in  computational  time.  For  example,  the  time  re¬ 
quired  to  solve  the  algebraic  equations  by  elimination  (at  a  fixed  x 
location)  is  proportional  to  the  number  of  unknowns,  whereas  the  time 
required  to  solve  the  equations  by  matrix  inversion  is  proportional  to 
the  square  of  the  number  of  unknowns  (Ref,  13).  On  the  IBM  370/155, 
a  solution  at  a  fixed  x  location  required  approximately  0.  6  sec  using 
200  grid  points  across  the  boundary  layer  (u  direction)  in  single  pre¬ 
cision.  The  corresponding  time  required  on  the  SDS  9300  was  approxi¬ 
mately  10  sec.  The  IBM  370/155  carried  8  digits  in  single  precision, 
whereas  the  SDS  9300  carried  12. 

The  general  method  of  solution  consisted  of  matching  the  inviscid 
and  viscous  flow  regions  in  the  nozzle  by  iterating  on  the  axial  pressure 
distribution.  An  initial  guess  of  the  axial  pressure  distribution  through¬ 
out  the  nozzle  is  made,  and  a  solution  is  obtained.  The  displacement 
thickness  calculated  in  this  solution  is  used  to  obtain  an  effective  nozzle 
geometry  which  in  turn  is  used  to  calculate  a  new  pressure  distribution 
from  one -dimensional,  perfect  gas,  expansion  theory.  A  typical  itera¬ 
tion  and  convergence  process  is  illustrated  in  Fig.  2.  The  use  of  one- 
dimensional  perfect  gas  expansion  theory  as  opposed  to  a  method-of- 
characteristics  solution,  for  example,  seems  justified  on  the  basis  of 
the  agreement  with  experimental  data  in  Section  V.  The  same  iteration 
process  as  used  here  and  some  suggestions  for  choosing  the  initial  pres¬ 
sure  distribution  are  discussed  further  by  Whitfield  and  Lewis  (Ref.  5). 

The  symbols  and  subroutines  used  in  the  Patankar  and  Spalding 
code  are  clearly  defined  and  discussed  in  Ref.  11.  Therefore,  the 
remainder  of  this  section  is  directed  toward  the  modifications  and 
additions  which  have  been  made  to  the  code  of  Ref.  11. 
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3.2  INPUT  CONDITIONS 


The  input  requirements  to  this  program  are  particularly  simple. 

The  input  was  modified  somewhat  from  that  of  the  original  code  (Ref.  11), 
and  the  input  variables  are  described  below  in  the  order  they  are  read- 
in  in  the  present  code  (Appendix  III), 

SYMBOL  DESCRIPTION 


KRAD 


IDIMEN 


NEQ 


KEX 


This  input  permits  the  treatment  of  plane  flows 
and  axisymmetric  flows  with  first-order  trans¬ 
verse  curvature.  Also,  although  not  pointed  out 
in  Ref.  11,  axisymmetric  flows  without  trans¬ 
verse  curvature  can  be  treated  by  setting 
KRAD  =  0.  Plane  flows  can  be  treated  by  set¬ 
ting  KRAD  =  0  and  rw(x)  =  1  (or  a  constant). 
Axisymmetric  flows  including  first-order  trans¬ 
verse  curvature  are  treated  by  setting  KRAD  =  1 
(at  least  not  zero)  and  using  the  actual  geometry 
rw(x). 

If  the  nozzle  considered  is  two-dimensional,  set 
IDIMEN  =  0.  This  sets  rw(x)  s  1  in  subroutine 
RAD.  If  the  nozzle  considered  is  axisymmetric 
set  IDIMEN  =  1. 

This  is  the  number  of  partial  differential  equations 
to  be  solved.  The  code  of  Ref.  11  also  includes 
the  solution  to  the  equation  for  the  conservation 
of  chemical  species.  However,  this  equation  was 
never  considered  in  this  investigation  and  only  the 
momentum  and  total  energy  equations  were  used, 
in  which  case  NEQ  =  2.  The  chemical  species 
equation  is,  however,  retained  in  the  present 
code. 

This  input  specifies  the  type  of  E  boundary.  It 
can  be  either  1,  2,  or  3,  according  to  whether 
the  E  boundary  is  a  wall,  free  boundary,  or  a 
symmetry  line,  respectively.  However,  in  the 
present  investigation  KEX  was  always  2,  i.  e. , 

E  was  a  free  boundary,  and  certain  modifications 
must  be  made  to  the  present  code  if  anything 
other  than  KEX  =  2  is  used. 
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KIN 


IHEAT 


N 


RE0RS 

ZETA 

PR(1) 

GAM 

ALPHA 

XR 

XL 

USUP 


This  is  similar  to  KEX  except  KIN  specified  the 
type  of  I  boundary.  For  the  present  code  KIN  =  1 
must  be  used;  otherwise  certain  modifications 
must  be  made. 

This  is  used  in  subroutine  FBC  and  indicates  a 
wall  temperature  boundary  condition  if  IHEAT  =  1 
and  a  wall  heat-transfer  rate  if  IHEAT  #  1.  In  the 
present  code,  the  only  wall  heat  flux  considered 
was  zero,  i.  e. ,  an  adiabatic  wall.  If  a  heat  flux 
other  than  zero  is  prescribed,  then  a  few  state¬ 
ments  in  subroutine  SLIP  must  be  modified  (see 
Ref.  11).  The  only  wall  temperatute  distribution 
considered  in  the  present  work  was  Tw/T0  =  con¬ 
stant.  However,  subroutine  FBC  is  easily  modi¬ 
fied  to  accommodate  any  desired  variation  of  Tw. 

This  is  the  number  of  strips  across  the  boundary 
layer,  i.  e. ,  in  the  to  direction.  It  must  always 
be  at  least  three  less  than  the  dimensions  of  the 
arrays  of  the  variables  across  the  layer,  e.  g. , 
the  maximum  N  which  can  be  used  with  the  di¬ 
mensions  of  the  arrays  in  the  present  code  is  197. 

Reservoir  Reynolds  number,  Re0jr*. 

Exponent  in  the  power  law  variation  of  viscosity 
with  temperature,  P 

Prandtl  number,  Pr 

Ratio  of  specific  heats,  7 

Nozzle  wall  half- angle,  a 

Longitudinal  radius  of_curvature  of  the  nozzle  up¬ 
stream  of  the  throat,  R 

Termination  condition  for  the  computations,  maxi¬ 
mum  value  of  x 

This  input  controls  the  location  of  the  E  boundary. 
It  is  associated  with  the  entrainment  rate  and  it 
will  be  discussed  in  the  following  subsection. 

USUP  was  varied  from  0.  99  to  0.  999,  and  was 
usually  0.  995. 
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YSTART  Initial  velocity  and  total  enthalpy  profiles  were 

calculated  from  the  expressions  u/uj^  =  2 and 
H/H0  =  Hw/H0  +  [  1-(HW/H0)  ] u/ u]£  where  n  =  y/yE 
and  y^  =  YSTART.  Suggestions  for  choosing 
YSTART  will  be  given  in  the  following  subsection. 
A  typical  value  of  YSTART  is  0.  5. 

TWT0  Ratio  of  wall  to  total  temperature,  TW/TQ 

XSTEP  Specification  of  the  integration  step-size  in  the 

x  direction  in  terms  of  local  wall  radius,  e.  g. , 
step  size,  DX,  is  given  by  the  product  of  XSTEP 
times  rw.  A  typical  value  of  XSTEP  is  0.  05. 


LMAX  Number  of  x  locations  where  an  input  pressure  is 

specified. 

XX(L)  Array  of  x  locations  where  an  input  pressure  is 

specified.  Array  goes  from  1  to  LMAX. 

P0P(L)  Input  pressure  for  corresponding  x  location,  XX(L). 

Array  goes  from  1  to  LMAX. 

The  last  three  inputs,  LMAX,  XX(L),  and  POP(L),  arc  read-in  in 
subroutine  PRE,  the  other  inputs  are  read -in  in  subroutine  BEGIN. 


Because  of  the  nature  of  the  process  involved  in  iterating  on  the 
axial  pressure  distribution,  it  is  important  to  input  a  smooth  initial 
pressure  distribution.  For  the  present  computations,  the  initial  pres¬ 
sure  distribution  was  calculated  using  the  actual  geometry  at  and  up¬ 
stream  of  the  throat  and  some  effective  inviscid  nozzle  wall  downstream 


of  the  throat.  (Actually,  the  slope  of  the  assumed  inviscid  nozzle  wall 
just  downstream  of  the  throat  was  matched  to  the  nozzle  wall  half-angle 
to  produce  a  smooth  wall  in  order  to  have  a  smooth  pressure  distribu¬ 
tion.  )  The  simple  program  used  to  calculate  the  initial  pressure  dis¬ 
tribution  for  the  present  computations  is  included  in  Appendix  IV  for 
convenience  of  users  where  such  an  approximation  is  adequate. 


The  particular  version  of  the  code  presented  in  Appendix  IV  uses 
6V  ~  x^/2  for  axisymmetrie  nozzles  and  S'1'  ~  x  for  two-dimensional 
nozzles.  Although  6:!<~  x^/2  was  not  used  to  calculate  the  initial 


pressure  distribution  for  all  axisymmetrie  nozzle  computations  pre¬ 


sented  herein,  it  appears  to  provide  a  reasonable  approximation  to 
the  variation  of  6*  in  axisymmetrie  nozzles. 
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Inasmuch  as  several  wind  tunnel  and  laser  nozzles  have  conical 
sections  downstream  of  the  throat  and  constant  longitudinal  radius  of 
curvature  for  the  upstream  converging  portions,  this  general  nozzle 
geometry  (which  can  be  sufficiently  described  by  the  inputs  ALPHA 
and  XE)  was  considered  in  this  investigation.  However,  subroutine 
RAD  can  be  easily  modified  to  include  any  geometry,  such  as  for 
example  contoured  nozzles  which  were  analyzed  in  Ref.  5  using  the 
same  governing  equations  as  used  here. 


3.3  SOME  POSSIBLE  COMPUTATIONAL  DIFFICULTIES 

The  mass -transfer  rate,  or  entrainment  rate,  across  the  E 
boundary  essentially  governs  the  location  of  the  edge  of  the  boundary 
layer.  This  technique  of  locating  the  edge  of  the  boundary  layer  has 
certain  advantages  in  analyzing  low -density  nozzle  flows;  however,  it 
might  also  cause  some  difficulties  if  not  handled  correctly.  The  en¬ 
trainment  rate  was  calculated  in  the  present  investigation  by  evaluating 
the  momentum  equation,  Eq.  (17),  along  a  constant  u  line,  denoted  as 
Ug,  near  the  E  boundary.  This  technique  is  discussed  in  Ref.  11.  The 
scheme  requires  the  specification  of  the  velocity  along  u  =  Ug  (where 
Ug  was  taken  as  0.  9)  at  the  next  downstream  station.  This  velocity  is 
denoted  as  ug  and  it  is  suggested  in  Ref.  11  that  it  be  taken  as 
ug  =  0.  99  ug,  where  ug  is  the  velocity  at  the  edge  of  the  boundary 
layer  at  the  next  downstream  station.  (Note  that  ug  can  be  calculated 
from  the  Euler  equation  since  dp/dx  is  presumed  known  for  the  par¬ 
ticular  iteration. )  It  was  found  in  the  present  work,  however,  that 
more  flexibility  could  be  obtained  with  the  program  if  ug  was  taken 
as  ug  =  (USUP)ug  and  USUP  was  input  for  each  solution.  The  quantity 
USUP  provides  a  means  of  suppressing  the  outer  edge  of  the  boundary 
layer  which  is  advantageous  in  treating  flows  which  are  nearly  merged. 
For  example,  during  the  process  of  iterating  on  the  axial  pressure  dis¬ 
tribution,  it  was  observed  (Fig.  2,  and  also  Ref.  5)  that  6*  resulting 
from  the  first  two  iterations  usually  provided  upper  and  lower  bounds 
on  the  final  converged  5'r'.  Therefore,  if  the  flow  is  sufficiently  rare¬ 
fied,  the  pressure  distribution  resulting  from  the  thinnest  6"  may  be 
such  that  the  following  iteration  would  predict  a  merged  flow  when  in 
fact  the  flow  is  not  merged  and  could  be  calculated  if  a  better  guess  for 
the  initial  pressure  distribution  could  be  made.  In  some  cases  the 
calculation  of  merged  flow  in  the  iteration  process  can  be  avoided  by 
using  a  small  value  of  USUP  (a  value  of  0.  99  is  herein  regarded  as  small 
and  0.  999  is  regarded  as  large)  to  suppress  the  edge  of  the  boundary 
layer  and  prevent  an  indication  of  merging.  The  suppression  seems  to 
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apply  only  the  outer  edge  of  the  layer  and  the  calculated  profiles  over 
most  of  the  layer  remain  essentially  unchanged.  For  this  reason,  the 
converged  6*,  or  pressure  distribution,  obtained  using  a  small  USUP 
usually  differs  by  a  small  amount  from  that  obtained  using  a  larger 
USUP,  say  0.  995.  After  convergence  is  obtained  using  a  small  value 
of  USUP  a  final  solution  can  be  obtained  using  a  larger  value.  The 
value  of  0.  99  for  ug/ug  as  suggested  in  Ref.  11  seems  to  excessively 
suppress  the  boundary  layer  for  the  present  internal  flow  calculations. 
For  most  solutions  reported  herein,  USUP  was  0.  995. 

The  step-size  along  the  x  component  was  taken  for  most  of  the 
solutions  as  5  or  7.  5  percent  of  the  local  wall  radius.  This  step-size 
was  sufficiently  small  for  most  problems.  However,  if  calculations 
of  properties  in  the  nozzle  throat  region  are  of  particular  interest, 
such  as  the  nozzle  wall  heat-transfer  rate,  a  smaller  step-size  may  be 
desirable. 

In  some  applications  where  nozzle  exit  properties  are  of  particu¬ 
lar  interest,  it  may  be  desirable  to  conserve  computer  time  and  use  a 
relatively  large  x  component  step-size.  For  such  problems  some  diffi¬ 
culties  might  be  encountered  in  starting  the  solutions.  Consider  a 
finite -difference  form  of  Eqs.  (19)  and  (20) 

^B-^D  =  ^E-^u  "  V  (24) 

If  Re0)  r*  and/or  YSTART  is  small,  then  (^E^'I^U  calculated  from  the 
initial  profile  will  be  small.  Depending  on  the  initial  profiles  and  flow 
conditions,  mg  may  be  positive  for  the  first  few  stations  and  therefore 
(^E-^'I)d  could  be  less  than  (ip E~^l)u  for  these  first  few  stations  and 
might  even  become  negative.  To  circumvent  this  difficulty  one  could 
reduce  the  step-size  (xq-xu)*  However,  if  this  is  not  desirable  in 
view  of  computational  time  requirements,  another  approach  is  to  in¬ 
crease  YSTART  in  order  to  increase  (^E^I^U  at  the  first  station.  The 
first  few  station  solutions  would  consequently  not  be  as  accurate  as 
solutions  obtained  by  using  a  small  x  step-size.  Therefore,  although 
the  solutions  are  started  at  the  beginning  of  the  converging  portion  of 
the  nozzle  upstream  of  the  throat,  this  technique  of  increasing  YSTART 
should  be  used  with  caution  if  accurate  solutions  (particularly  proper¬ 
ties  which  depend  on  derivatives  near  the  wall  such  as  heat-transfer 
rate)  in  the  throat  region  are  required.  As  stated  previously,  a  typical 
value  of  YSTART  was  0.  5.  The  largest  value  of  YSTART  used  to  obtain 
solutions  was  unity,  but  this  depends  on  XR. 
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Heat  transfer  and  viscous  effects  change  the  effective  nozzle  throat 
from  the  location  corresponding  to  the  actual  geometric  nozzle  throat. 
The  subroutine  NEWPP0  searches  for  the  minimum  effective  nozzle 
radius  and  uses  it  for  the  throat  in  calculating  a  one -dimensional  pres¬ 
sure  distribution  for  the  next  iteration.  However,  if  a  sufficient  num¬ 
ber  of  solutions  near  the  throat  are  not  taken,  then  the  true  minimum 
area  and  its  location  used  in  calculating  a  new  pressure  distribution 
may  not  be  accurately  approximated.  Then,  the  resulting  dp/dx  in  the 
throat  region  for  the  following  iteration  may  not  be  smooth.  In  this 
case,  the  error  in  dp/dx  for  each  successive  iteration  would  become 
worse.  This  is  especially  a  problem  in  solutions  for  flows  which  have 
adiabatic  or  relatively  hot  nozzle  walls.  It  is  suggested  that  if  the 
initial  guess  of  the  distribution  is  not  a  particularly  good  one,  e.g. , 
if  the  calculated  6',:  at  the  exit  is  not  within  20  to  30  percent  of  the 
assumed  6'“',  then  an  improved  pressure  distribution  should  be  calcu¬ 
lated  by  the  program  in  Appendix  IV  (or  some  similar  method).  This 
ensures  a  smooth  pressure  distribution,  and  since  convergence  is,  in 
general,  much  faster  in  the  throat  region  than  further  downstream 
where  the  relative  displacement  thicknesses  are  larger,  few,  if  any, 
extra  iterations  are  required. 

Although  areas  have  been  pointed  out  where  computational  difficul¬ 
ties  have  been  encountered  with  this  program,  it  should  also  be  pointed 
out  that  this  numerical  scheme  is  actually  rather  rugged,  as  for 
example  compared  to  the  method  of  Ref.  5.  It  is  in  general  not  sensi¬ 
tive  to  input  conditions  and  seldom  "blows  up.  " 


SECTION  IV 
NUMERICAL  RESULTS 


Solutions  are  presented  in  this  section  for  typical  nozzle  configura¬ 
tions  which  provide  some  indication  of  the  effect  of  Re0  r*,  7,  Pr,  ?, 
Tw/T0,  a,  transverse  curvature,  and  two-dimensional  versus  axisym- 
metric  flows.  Table  I  (Appendix  II)  summarizes  the  conditions  of  the 
solutions  presented  in  this  section. 

Of  particular  interest  in  nozzle  flows  is  the  displacement  thickness, 
S'".  By  using  6*  and  the  nozzle  geometry,  an  effective  inviscid  nozzle 
radius  can  be  determined  from  which  Mach  number  and  other  flow  prop¬ 
erties  outside  the  boundary  layer  can  be  estimated.  The  displacement 
thicknesses  for  Conditions  1  to  8  of  Table  I  are  presented  in  Figs.  3  to  10. 
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Solutions  were  started  at  the  beginning  of  the  converging  portion  of  the 
nozzles  where  z/r!,!  =  -3.  However,  to  conserve  computing  time,  ex¬ 
tremely  small  step-sizes  were  not  utilized  for  all  the  solutions  in  this 
region  (as  discussed  in  Section  III),  and  results  are  not  presented  in 
Figs.  3  to  10  for  z/r'1'  <  -2.  Some  solutions  were  repeated,  however, 
with  smaller  step-sizes  with  no  appreciable  changes  in  the  results  pre¬ 
sented. 

At  least  for  some  values  of  ReQ  r*  in  Figs.  3  to  7  and  9  to  10, 
negative  displacement  thicknesses  were  calculated  in  the  throat  region. 
This  is  due  to  the  relatively  cool  wall  increasing  the  gas  density  and 
hence  the  local  mass  flux  near  the  wall.  The  adiabatic  wall  results 
presented  in  Fig.  8  do  not  indicate  negative  6''  for  the  same  flow  con¬ 
ditions.  Similar  results  concerning  the  calculation  of  negative  6;,:  have 
been  reported  previously,  e.  g. ,  Potter  and  Carden  (Ref.  3)  and 
Whitfield  and  Lewis  (Refs.  4  and  5). 

Some  indication  of  the  effect  of  using  various  gases  in  a  fixed  nozzle 
geometry  with  fixed  flow  conditions  is  provided  in  Figs.  3  to  5.  The 
specific  heat  ratio,  Prandtl  number,  and  power  law  viscosity  variation 
with  temperature  of  Figs.  3  to  5  closely  approximate  the  properties  of 
carbon  dioxide  (C02),  helium  (He),  and  nitrogen  (N2),  respectively. 

For  example,  by  considering  the  nozzle  exit  displacement  thickness  for 
Re0j  r*  =  104.  one  observes  that,  for  the  indicated  flow  conditions,  6;-: 
using  He  is  about  200  percent  of  that  when  using  C02,  and  6'"'  using  N2  is 
about  150  percent  of  that  when  using  C02. 

Results  are  presented  in  Figs.  5  and  6  for  identical  conditions 
except  for  the  nozzle  wall  half-angle.  Figure  5  has  a  wall  half-angle 
of  10-deg,  and  Fig.  0  has  a  half -angle  of  15  deg.  For  the  same  geo¬ 
metric  area  ratio,  the  expansion  process  is  more  rapid  for  the  a  =  15 
deg  nozzle  than  the  a  =  10  deg.  Also,  the  nozzle  wall  length  for  the 
same  area  ratio  is  longer  for  the  10-deg  nozzle  than  for  the  15 -deg 
nozzle.  The  result  is  to  produce  a  larger  exit  6"  for  the  10-deg  nozzle 
than  for  the  15-deg  nozzle  for  the  same  geometric  area  ratio.  For 
Re0  r*  =  3  x  103,  the  10-deg  nozzle  exit  6  ’  is  about  18  percent  larger 
than  that  of  the  15-deg  nozzle. 

It  should  be  pointed  out  that  although  6'1'  is  relatively  small  because 
of  the  cool  walls,  5  is  not  necessarily  small.  For  example,  in  Fig.  5 
for  ReQ  *  =  10^,  the  flow  merges  (i.  e.  ,  the  boundary  layer  completely 
fills  the*  nozzle)  at  a  point  where  6/r*  ~3.  For  these  conditions, 

6'‘'/6  ~ 0(1/  10).  The  ratio  6 '■' / 6  tends  to  increase  with  wall  temperature. 
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In  general,  the  wall  temperature  seems  to  have  a  stronger  effect  on  6'"' 
than  on  6  (Ref.  4).  In  Figs.  5,  7,  and  8  the  effect  of  nozzle  wall  tem¬ 
perature  was  investigated,  with  other  conditions  held  constant,  by  con¬ 
sidering  Tw/T0  =  1/10,  Tw/T0  =  1/3,  and  an  adiabatic  wall.  For 
Re0,  r*  =  10^,  the  ratio  6*/ 6  at  the  nozzle  exit  was  0,30,  0.48,  and 
0.  68,  respectively. 

The  effect  of  Prandtl  number  was  investigated  by  repeating  the  con¬ 
ditions  of  Fig.  5  but  with  Pr  =  1  in  place  of  5/7.  The  results  are  pre¬ 
sented  in  Fig.  9.  The  displacement  thickness  was  found  to  increase 
for  Pr  =  1  by  approximately  20  percent  at  the  nozzle  exit. 

The  conditions  of  Fig.  5  were  also  repeated  using  £  =  1  in  place 
of  £  =  2/3  to  investigate  the  effect  of  the  power  law  variation  of  vicosity 
with  temperature.  For  £  =  1  the  displacement  thickness  was  found  to 
decrease  by  20  to  25  percent  below  that  for  ^  =2/3  (see  Figs.  5  and  10). 

_ f  _ 

Note  that,  since  m  =  T  *and  T  <  1,  then  the  flow  is  more  viscous  for 
£  =  2/3  than  for  £  =  1. 

Velocity  and  temperature  profiles  calculated  with  and  without  trans¬ 
verse  curvature  terms  are  presented  in  Fig.  11  for  Conditions  3  and  9 
(Table  I)  with  ReQ  r*  =  3  x  103.  Neglecting  transverse  curvature  de¬ 
creases  the  boundary -layer  thickness  for  internal  flow.  The  effect  of 
transverse  curvature  is  negligible  for  thin  boundary  layers,  but  may  be 
significant  for  thick  boundary  layers.  Further  results  and  discussion 
concerning  the  effects  of  transverse  curvature  are  reported  in  Ref.  5. 

The  displacement  thicknesses  of  an  axisymmetric  and  two- 
dimensional  nozzle  are  presented  in  Fig.  12.  For  each  nozzle, 

Ae/A*  =  5,  and  therefore,  the  two-dimensional  nozzle  is  considerably 
longer  since  a  =  10  deg  for  each.  The  result  is  to  produce  a  signifi¬ 
cantly  larger  displacement  thickness  for  the  two-dimensional  nozzle 
than  for  the  axisymmetric  nozzle. 


SECTION  V 

COMPARISONS  WITH  EXPERIMENTAL  DATA  AND 
OTHER  THEORETICAL  RESULTS 


Curves  of  Re0,  r*/(Por>i!)  versus  total  temperature  for  various  gases 
are  presented  in  Fig.  13  for  the  convenience  of  determining  Re0<  r*.  The 
viscosities  used  for  all  gases  in  Fig.  13  except  carbon  dioxide  (CO2)  are 
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those  given  by  Svehla  (Ref.  14).  The  viscosity  of  CO2  was  taken  from 
Table  8.4-2  of  Hirschfelder,  Curtiss,  and  Bird  (Ref.  15). 

Highly  viscous  nozzle  flows  are  usually  associated  with  nozzles  of 
small  physical  size.  The  spatial  resolution  of  experimental  measure¬ 
ments  in  such  nozzles  is  obviously  restricted.  However,  highly  viscous 
flows  can  also  be  produced  in  large  nozzles  if  sufficiently  large  volume 
flows  can  be  pumped  at  low  pressures,  thereby  permitting  more  detailed 
investigations  of  the  boundary  layer.  Such  investigations  are  possible 
in  the  Aerospace  Research  Chamber  10V  in  the  VKF  at  AEDC.  The 
present  and  some  previous  calculation  methods  will  be  compared  with 
experimental  data  taken  in  Chamber  10V  using  a  nominal  Mach  three 
nozzle,  denoted  M3  nozzle.  This  is  a  10 -deg  half-angle  conical  nozzle 
with  d*  =27  cm,  de  =  76.  2  cm,  and  R  =  11.2  cm.  The  walls  of  this 
nozzle  were  cooled  with  liquid  nitrogen  to  maintain  a  constant  nozzle 
wall  temperature  of  about  85°K. 

Calculated  and  measured  pitot  pressure  profiles  at  the  nozzle  throat 
are  presented  in  Fig.  14.  Nitrogen  was  the  test  gas  used  for  these  and 
all  other  experimental  data  presented  in  this  report.  The  present  results 
are  in  good  agreement  with  the  results  using  the  method  of  Ref.  5.  The 
calculated  profiles  are  in  relatively  good  agreement  with  the  measured 
profile  for  these  conditions. 

Calculated  and  measured  pitot  pressure  profiles  at  the  nozzle  exit 
are  presented  in  Fig.  15.  This  calculation  was  performed  using  £  =  2/3, 
whereas  some  other  results  presented  herein,  e.  g. ,  Fig.  14,  were  ob¬ 
tained  using  Sutherland's  viscosity  law  for  nitrogen.  However,  over  the 
temperature  range  of  the  experimental  data  taken  in  Chamber  10V,  there 
was  negligible  difference  in  the  results  using  either  viscosity  variation 
with  temperature. 

Present  calculations  are  compared  to  calculations  by  the  method  of 
Rae  (Ref.  10)  and  experimental  data  in  Figs.  16  and  17.  The  present 
results  are  in  relatively  good  agreement  with  the  experimental  data, 
but  Rae's  method  is  found  to  overestimate  the  size  of  the  viscous  region 
for  these  conditions.  Results  were  also  presented  in  Ref.  5  for  the 
conditions  of  Fig.  17.  The  calculated  profile  from  Ref.  5  is  almost 
identical  with  the  present  result  in  Fig.  17  (see  Ref.  5)  and  is  not  pre¬ 
sented. 

Pitot  pressure  profiles  at  about  7  throat  radii  downstream  of  the 
M3  nozzle,  throat  are  presented  in  Fig.  18  for  Re0  r*  =  1170.  A  rela¬ 
tive  minimum  exists  at  the  nozzle  centerline  in  the  experimentally 
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measured  profile,  and  slight  "humps"  or  relative  maximums  exist  near 
the  edge  of  the  boundary  layer.  It  should  be  pointed  out  that  the  exist¬ 
ence  of  humps  does  not  necessarily  imply  that  the  flow  is  non- 
isentropic.  Imagine  a  radial  pitot  pressure  profile  in  an  inviscid  con¬ 
toured  supersonic  nozzle  which  is  sufficiently  far  downstream  of  the 
throat  to  pass  through  a  portion  of  the  uniform  parallel  flow  near  the 
centerline.  Such  a  pitot  profile  would  have  lower  values  near  the 
nozzle  centerline  than  near  the  wall  because  of  the  larger  Mach  num¬ 
bers  near  the  centerline.  However,  the  matching  of  such  an  inviscid 
flow  with  a  realistic  viscous  flow  requires  the  pitot  pressure  near  the 
wall  to  approach  the  local  static  pressure  which  is  less  than  the  center- 
line  pitot  pressure.  Therefore,  it  is  not  difficult  to  conceive  of  humps 
in  such  a  radial  pitot  profile  because  of  the  matching  of  inviscid  and 
viscous  flow.  This  argument  is  based  on  flow  in  a  contoured  nozzle. 

It  is  applied  to  the  present  case  because  the  displacement  thickness 
effectively  contours  the  nozzle.  A  more  accurate  approach  of  investi¬ 
gating  such  flows  would  be  to  remove  the  present  assumption  of  one- 
dimensional  inviscid  flow  outside  the  boundary  layer  and  obtain  more 
accurate  solutions,  such  as,  inviscid  method -of -characteristics  solu¬ 
tions.  However,  the  pitot  pressure  profile  is  relatively  well  predicted 
in  Fig.  18,  and  for  this  work,  the  assumption  of  one -dimensional  in¬ 
viscid  flow  outside  the  boundary  layer  is  considered  acceptable. 

Just  as  radial  pitot  profiles  with  humps  do  not  necessarily  imply 
that  the  flow  is  non-isentropic,  aflat  radial  profile  does  not  necessarily 
imply  that  the  flow  is  isentropic.  Consider  the  pitot  pressure  data  in 
Fig.  19.  Both  the  calculated  and  measured  profiles  are  relatively  flat 
for  y/rw  larger  than  about  0.  7.  However,  in  this  case  the  pitot  pro¬ 
file  is  not  a  good  measure  of  the  extent  of  the  nozzle  wall  viscous 
effects  as  shown  in  Fig.  20.  From  Fig.  20,  the  bound  ary -layer  thick¬ 
ness  as  estimated  from  the  pitot  profile  is  about  55  percent  of  the  nozzle 
radius,  whereas  it  is  calculated  to  be  actually  over  80  percent.  The 
pitot  profile  for  the  case  in  Fig.  20  implies  that  20  percent  of  the  cross- 
sectional  area  of  the  nozzle  at  this  point  is  core  flow,  whereas  actually 
less  than  4  percent  is  core  flow.  The  reason  a  pitot  profile  might  lead 
one  astray  is  associated  with  the  temperature  or  thermal  boundary 
layer  since  pitot  pressure  depends,  among  other  things,  onu/(T)*'2. 
The  velocity  variation  is  usually  well  behaved  and  fairly  accurately 
predicted  by  simple  analytical  expressions  based  on  boundary-layer 
thickness  (Ref.  4);  however,  this  is  not  the  case  with  the  thermal 
boundary  layer.  The  temperature  variation  depends  not  only  on  local 
wall  and  edge  values  and  gradients  but  also  on  the  upstream  conditions. 
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Therefore,  some  consideration  should  be  given  to  the  temperature  vari¬ 
ation  in  order  to  place  limits -of -confidence  on  the  use  of  pitot  pressure 
as  an  indication  of  the  nozzle  wall  viscous  effects. 

Carden  (Ref.  16)  measured  local  heat-transfer  coefficients  in  an 
axisymmetric  nozzle  for  Re0<ri!J  =  5  x  102.  The  experimental  data 
were  compared  with  calculated  heat -transfer  coefficients  in  Refs.  5 
and  16.  However,  Whitfield  made  a  mistake  in  Ref.  5  and  presented 
solutions  for  r'1'  =  0.262  cm  instead  of  d'*5  =  0.  262  cm  (0.  103  in. )  which 
corresponds  to  the  actual  nozzle  throat  dimension.  Because  of  this 
error,  the  calculated  heat-transfer  coefficients  (Ref.  5)  from  the 
iterated  solutions  were  significantly  below  the  experimental  data  of 
Carden  (Ref.  16).  Heat -transfer  coefficients  were  calculated  using 
the  present  method  but  with  the  inappropriate  throat  dimension  of 
r*  =  0.  262  cm,  and  good  agreement  was  obtained  with  the  calculated 
results  of  Ref.  5.  Results  were  also  obtained  with  the  present  method 
using  the  proper  throat  dimension  of  d5,:  =  0.  262  cm,  and  good  agree¬ 
ment  with  the  experimental  data  was  obtained  as  shown  in  Fig.  21. 

Also,  in  Fig.  21  are  results  from  two  of  the  methods  Carden  (Ref.  16) 
used  for  calculating  the  heat -transfer  coefficient  and  one  solution  from 
Ref.  5.  Although  the  result  from  Ref.  5  which  is  presented  in  Fig.  21 
was  not  iterated  to  include  the  higher-order  displacement  effect,  the 
result  was  obtained  using  the  pressure  distribution  corresponding  to  the 
proper  nozzle  geometry  and,  therefore,  is  included  in  Fig.  21.  This 
result  from  Ref.  5  is  in  relatively  good  agreement  with  the  present 
result.  All  calculation  methods  underestimate  the  most  upstream  un¬ 
corrected  experimental  data  point  in  Fig.  21.  However,  Carden  points 
out  that  radiation  from  the  arc  used  to  heat  the  gas  could  increase  the 
total  heat-transfer  rate  to  this  portion  of  the  nozzle.  No  corrections  for 
radiation  heat  transfer  were  made  to  the  experimental  data. 

It  might  be  pointed  out  that  the  heat-transfer  rate  printed  out  in  the 
present  program,  denoted  as  AJI(l),  is  equal  to  qw/(p0H0^^2). 

Although  made  dimensionless,  the  numerical  scheme  used  in  Ref.  11  to 
calculate  AJI(l)  is  retained  in  the  present  code  in  subroutine  WALL. 

SECTION  VI 
CONCLUSIONS 

The  present  results  were  shown  to  be  in  good  agreement  with  re¬ 
sults  from  the  method  of  Ref.  5,  Although  it  was  noted  during  the  pres¬ 
ent  investigation  that  more  axial  stations  were  required  using  the 
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present  program  than  that  of  Ref.  5  to  obtain  the  same  degree  of  accu¬ 
racy  (presumably  this  is  due  to  the  transformation  variables  used  in 
Ref.  5,  which  are  advantageous  for  thin  boundary  layers  but  disadvan¬ 
tageous  for  thick  internal  boundary  layers),  the  present  program  is 
much  smaller,  faster  and  easier  to  use.  Comparisons  with  results 
from  Rae's  method  (Ref.  10)  and  experimental  data  indicated  the  pres¬ 
ent  method  to  be  more  accurate  than  Rae's  for  the  conditions  considered. 
It  should  be  remembered,  however,  that  the  present  work  is  not  applic¬ 
able  to  merged  flows  without  modification,  whereas  Rae's  method  is 
more  applicable  to  merged  flows. 

Consistent  agreement  was  obtained  between  the  present  results 
and  experimental  data.  Except  for  the  heat-transfer  data,  the  experi¬ 
mental  data  were  taken  in  a  relatively  large  nozzle  where  the  boundary 
layer  could  be  accurately  probed.  Similar  agreement  was  obtained  in 
Ref.  2  where  extensive  data  were  obtained  in  small  nozzles, 
r*  ~0  (0.  1  cm).  The  numerical  method  used  in  Ref.  2  was  developed 
by  Mayne  (Ref.  2),  as  stated  previously,  and  differs  from  the  original 
Patankar  and  Spalding  method  (Ref.  11)  primarily  in  the  numerical 
scheme  near  the  wall. 

The  present  method  has  not  been  modified  for  use  in  nozzle  design; 
however,  it  is  a  straightforward  matter  to  make  such  a  modification 
(see  Ref.  5).  To  design  a  nozzle,  a  new  nozzle  geometry,  rw(x),  is 
calculated  from  the  displacement  thickness  from  each  iteration,  and 
the  new  geometry  is  used  for  the  next  iteration  while  the  desired  axial 
pressure  distribution  is  maintained  for  each  iteration.  It  was  found  in 
Ref.  5  that  convergence  was,  in  general,  more  rapid  in  iterating  on 
the  nozzle  geometry  than  on  the  axial  pressure  distribution. 
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Fig.  1  Definition  of  Coordinate  System 
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Fig.  2  Displacement  Thickness  of  Successive  Iterations 
for  Condition  No.  5  (Table  I) 
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Fig.  4  Displacement  Thickness  for  Condition  No.  2  (Table  1} 
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Fig.  5  Displacement  Thickness  for  Condition  No.  3  (Table  I) 
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Fig.  6  Displacement  Thickness  for  Condition  No.  4  (Table  1} 
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Fig.  8  Displacement  Thickness  for  Condition  No.  6  (Table  I) 
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Fig.  9  Displacement  Thickness  for  Condition  No.  7  (Table  I) 
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Fig.  10  Displacement  Thickness  for  Condition  No.  8  (Table  I) 
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Condition  No.  3  (Table  I), 
with  Transverse  Curvature 


Fig.  11  Nozzle  Exit  Velocity  and  Temperature  Profiles 
for  Condition  Nos.  3  and  9  (Table  I) 
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Fig.  12  Displacement  Thicknesses  of  an  Axisymmetric  and 
Two-Dimensional  Nozzle  with  Equal  Area  Ratios 
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Fig.  13  Reservoir  Reynolds  Number  as  a  Function 
of  Pb.  T0,  and  r*  for  Various  Gases 
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Fig.  14  Calculated  and  Measured  Pitot 
at  the  M3  Nozzle  Throat  for  F 
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Fig.  16  Calculated  and  Measured  Pitot  Pressure  Profiles 
at  the  M3  Nozzle  Exit  for  Ret,,r«  =  900 


42 


AEDC-TR-73-52 


Sym 

PQ,  torr 

T  °K 

o» 

T  ,  °K 
w  * 

Source 

O 

0.50 

280 

89 

Experimental  Data 

0.50 

300 

100 

Calculated 

— 

0.50 

300 

89 

Calculated 

0  0.2  0.4  0.6  0.8  1.0 


y/re 

Fig.  17  Calculated  and  Measured  Pitot  Pressure  Profiles 
at  the  M3  Nozzle  Exit  for  Re0#r*  =  4500 
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Fig.  18  Calculated  and  Measured  Pitot  Pressure  Profiles  Downstream 
of  the  M3  Nozzle  Throat  for  Reo.r*  “  1170 
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Fig.  19  Calculated  and  Measured  Pitot  Pressure  Profiles  Downstream 
of  the  M3  Nozzle  Throat  for  Re^,r»  =  330 
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Fig.  20  Calculated  Velocity,  Temperature,  and  Pitot  Pressure  Profiles 
Downstream  of  the  M3  Nozzle  Throat  for  Re,,  r.  =  330 
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Fig.  21  Calculated  and  Measured  Heat-Transfer  Coefficients 
in  a  Low-Density  Nozzle 
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TABLE  I 

CONDITIONS  FOR  THE  NUMERICAL  SOLUTIONS  PRESENTED  IN  SECTION  IV 


Condition 

No. 

D 

Pr 

n 

Tw 

/T0 

a, 

deg 

A/ A* 

R 

Reo,  r* 

Comment 

1 

9/7 

3/4 

3/4 

1/10 

10 

25 

3 

3  x  103,  104,  105 

Axisymmetri 

c  Nozzle  with 

Transverse  Curvature 

2 

5/3 

2/3 

5/7 

10 

5  x  103,  104,  2  x  104,  5  x  104.  105 

3 

7/5 

5/7 

2/3 

10 

103,  3  x  103,  104,  3  x  104,  105 

4 

■B 

3  x  103,  104,  10s 

5 

1/3 

99 

5  x  103,  104,  3  x  104,  105 

6 

Adiabatic 

104,  105 

7 

1 

1/10 

3  x  103,  104,  3  x  104,  105 

8 

5/7 

1 

103,  104,  105 

9 

5/7 

2/3 

■ 

I 

3  x  103 

Axisymmetri 

c  Nozzle  with- 

I 

out  Transverse  Curvature 

10 

5/7 

2/3 

■ 

1 

3  x  103 

Two-Dimensional  Nozzle 
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APPENDIX  III 

BOUNDARY-LAYER  COMPUTER  CODE 


FORTRAN  IV  G  LEVEL  20 


MAIN 


DATE  ■  72286 


22/52/18 


C  MAIN  A  2 

OOOl  COMMON  /GEN/  PEI ,AMI ,AHE, 0PDX.PREF(2I ,PR(2» >P(2I , DEN, AMUt XU.XD.XP.  A  3 

. 1XL »  OX*  INI  GiCSALFaValPHA  |XR  (RE’bRSiGAMt’ZET  A,PPd,TWTO,YSTAHT,USUP«IOI  **A  V 

_ 2MEN,lHEATt2 tT0,XSTEP/l/N,NPl,NP2,NP3,NEQiNPHtKEX,KIN,KASE,KRAO/ B/B  A  5 


3ETA ,  GAMA ( 2  1 1 T AUI ,  T AUE i AJI (21 , AJEI2I , 1N0I (2) . I  N0EC2I /V/UI200 1 , F< 2, 2  A 
400 ) ,R ( 2001 ,RHO(2001 ,0M( 200) , Y(200) /C /SC ( 200 > , AU< 200 > , BU (200 1 , CU( 20  A 
50 )  i A(  2>200VtB (2't200T>  ,C (2,'260l7d/YR('200 lVUR(200 1,RR(200)  ,HR«  200 1  iXH'A 


-  0002 

— 

6(200) # P I TOT ( 200 ) *  TEMP ( 200 ) /E /DSTAR (300) »  XRS( 300) (RWRS (300) t COSAL (3 
700) 

CDMHON  /L/  AK.ALMG 

A 

..... 

A 

9 

“id* 

u 

0003 

i 

CONTINUE 

A 

12 

0004 

AK*1« 

A 

13 

0005 

INTG*0 

A 

14 

0006 

CALL  BEGIN 

A 

IS 

'  0007 

CALL  PRE  ( XU.DPOX'I 

A 

16 

0003 

XTH0AT»3. 141S9*XR/2. 

A 

17 

0009 

AM 1*0. 

nr 

IB 

0010 

AME-O. 

A 

19 

0011 

GO  TO  3 

A 

20 

0012 

2 

CALL  READY 

A 

21 

C 

THE  FOLLOWING  STATEMENT  IS  TO  KILL  SHOT  IF  PROBLEMS  (LIKE  MERGING) 

A 

22 

c 

OCCUR  IN  SUUROUTI NE  READY 

A 

23 

0013 

IF  (XO.GT.XL)  GO  TO  13 

A 

24 

0014 

. 3. 

CONTINUE 

A 

25 

0015 

•  INTG* INTG* 1 

A 

26 

C 

CHOICE  OF  FORWARO  STEP 

A 

27 

0014 

DX«XSTEP*R( 1) 

A 

2B 

0017 

XO-XU+DX 

29 

0019 

IF  ( .NOT. (XU. LT. XTHOAT •  ANO. XO.GT. XTHOAT ) )  GO  TO  4 

A 

30 

0019 

XD-XTHOAT 

A 

31 

0020 

DX*XD-XU 

A 

32 

0021 

J 

ITHOA  T“INTG*1 

A 

33 

0022 

4 

CONTINUE 

A 

34 

0023 

A 

35 

0024 

XO-XL 

A 

36 

0025 

OX-XD-XU 

A 

.37.,. 

002b 

5 

CONTINUE 

A 

3B 

0027 

CALL  PRE  (XOtOPOX) 

A 

39 

0028 

CALL  ENTRN 

A 

40 

0029 

CALL  PRE  (XU.OPOX) 

A 

41 

0030 

IF  (KASE.EQ.2)  GO  TO  6 

A 

42 

0031 

IF  (K1N.GQ.IJ  CALL  MASS  (XU,XO,AMI) 

A 

43 

0032 

IF  (KEX.EO. 1)  CALL  HASS  (XU.XO.AME) 

A 

44 

0033 

CALL  WALL 

A 

45 

0034 

6 

CALL  OUTPUT 

A 

46 

0035 

CALL  PRE  (XOiOPOX) 

A 

47 

0036 

CALL  COEFF 

nr 

48 

C 

SETTING  UP  VELOCITIES  AT  A  FREE  BOUNDARY 

A 

49 

C 

MODIFIED  FOLLOWING  STATEMENT  FOR  INTERNAL  CORf  FLOW 

A 

50 

0037 

IF  (KEX.EO. 2)  U(NP3)-SQRT(2.l*SURT(l.-PP0**<(GAM~l.J/GAMn 

A 

51 

0038 

"  IF '  (KIN.  E0.  2 )  U ( 1 1  “SORT (U ( 1 )  *U(  l  )-2T*(X6-XU)*0'P0X'/RM6T £TI 

A 

"52“' 

0039 

CALL  SOLVE  ( AU.BUtCUi U.NP3) 

A 

53 

c 

SETTING  UP  VELOCITIES  At  A  SYMMflRY  LINE 

nr 

sv~ 

0040 

IF  (KIN.NE.3)  GO  TO  7 

A 

55 

U( 1)*U( 2)  . 

A 

56 

0042 


IF  (KRAD.EO.O)  U(1)-.75*U(2)*.25*U(3I 


A  57 
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rimRAN  iv  c 

LEVEL 

20  MAIN  DATE  «  72286  22/52/18 

0043 

7 

IF  (KEX.E0.3I  U<  NP3 1«.  75*U(NP*2) ♦. 25*U(NP1 1 

"58 

0044 

IF  (NEO.EO. 11  GO  TO  14 

A 

$9 

004$ 

00  13  J-l.NPH 

A 

60 

0046 

00  8  I«2»NP2 

A 

61 

0047 

AU(II-A(J,II 

62 

004B 

BUm-BI  Jtl  J 

A 

63 

004° 

a 

cum-cu.n 

A 

64 

0050 

00  9  I-1.NP3 

A 

6$ 

0051 

9 

scin-F(jtn 

A 

66 

.  0052 

• 

CALL  SOLVE  I AU, BU.CU. SC ,NP31 

A 

67 

0053 

00  10  l>l,NP3 

A 

68 

0054 

10 

Fu.n-scm 

A 

69 

005$ 

IF  (KASE.E0.2)  GO  TO  11 

A 

70 

c 

SETTING  UP  HALL  VALUES  OF  F 

A 

71 

0036 

IF  (K1N.EQ.UAN0.1N0MJI.E0.2>  P'ij|l>-ni;VtteTA4GANAlJl|6>(j|2t-<r 

A 

"  72"'* 

l.*BETA-GAHA(J||*FU,3)l*.5/GANAIJI 

A 

73 

0057 

iF  (KEX.EO.l.AND. I NOE (J 1 . EU. 2 i  F ( j,NP3 )«( ( l.*BETA«-GANA( Jt l*F( J,NP2 

A 

~74 

11 -( 1 . ♦3ETA-GANA IJ II *F (J ,NPl)i *. 5/GAMAI J) 

A 

75 

c 

SETTING  UP  SYMMETRY-LINE  VALUES  OF  F 

A 

76  '* 

0058 

11 

IF  (KIN.NE.3I  GO  TO  12 

A 

77 

0050 

F(J,1I-F(J,2| 

A 

"78  ' 

0060 

IF  (KRAO.Ed.OI  FIJ.1I».75*F(J|2I+.25*F(J|3I 

A 

79 

0061 

12 

IF  (KEX.cQ.il  F ( J »NPJ 1*. 1 5*F ( Ji NP2)-».  25*F  ( J.NP1 ) 

A 

80 

0062 

13 

CONTINUE 

A 

81 

0063 

14 

XP-XU 

A 

82 

0064 

XPS( INTGI-XU 

A 

83 

0065 

RWPSC INTGI-flt 11 

A 

84 

0066 

COSAL( INTGI»CSALFA 

A 

85 

0067 

XU-XO 

A 

86 

0068 

PEI»PEI*0X*(R(1I*AMI-R(NP31*ANEI 

A 

87 

C 

THE  TERMINATION  CONDITION 

A 

88 

0069 

l . . 

IF  (XU.LT.XLI  GO  TO  2 

A 

89 

0070 

CALL  NEHPPO 

A 

90 

0071 

STOP 

A 

91 

0072 

END 

A 

92- 
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FOR THAN  IV  6  LEVEL  20  NEWPPO  DATE  *  722B6  22/52/18 


0001 

...  0002... 

SUBROUTINE  NEWPPO 

COMMON  /GEN/  PE  It AMI , AME, OPDA ,PREF 12 I ,PR (21 , P 1 21 ,DEN, AMU.XU.XO.XP, 
lxl. OX, INTG.CSALFA, ALPHA fXR.REORSf CAM, ZET  A,PP'0,THTQ(VSf  ART.USUP, IDI 
2M£N,IHEAT,Z,T0,XSTEP/£/DSTAR(300l .XRSI300I ,KWRS(300I ,C0SAL(300I 

B 

B 

"B" 

B 

1 

2 

'"3""' 

4 

“5  — 
6 

7 

8 

9 

10 

0003 

0004 

"  0005'" 
0006 
0007 
0008 

. 

0 IMF.NS ION  FMMI300I,  POPI300I,  RIRSI300) 

IF  I  1DIMEN.EQ.0I  RETURN 

G=GAM 

C 1= I G« 1.  1/2. 

C2-IG-1. 1/2. 

C 3=  1 0*1.  )  / 1  2 . *  1 G—  1 .  II 

B 

B 

"B  " 

8 

'  V 

B 

“ ' "oooh 

C4«C1**C3 

8 

11 

0010 

C5«C4*C2*C3 *2. 

B 

12 

0011 

DO  1  1-1, INTO 

B 

13 

0012 

1 

RIPSMI-RMRSUI-COSALIIXDSTARIII 

B 

14 

0013 

DO  2  1=1, INTG 

"B" 

15 

0014 

IF  (RIRSIII.LT. RIRSIMDI  GO  TO  3 

8 

16 

0015 

2 

CONTINUE 

8 

“17 

0016 

3 

RTHOAT-KMRSI Il-COSALI II*DSTARII 1 

8 

18 

0017 

DO  14  1=1 , INTG 

B 

19 

0018 

A»(RTHOAT /RIRSI I) 1**2 

8 

20 

0019 

IF  t  1-  1 1  4,4,6 

8 

21 

0020 

4 

FM1-0.5 

B 

22 

0021 

5 

B=( l.«C2*FMl*«2l 

6 

23 

0022 

FM-FH1-(A«B**(C3«1. I-C4*FM1*BI/IC5*FM1**2-C4*BI 

B 

24 

0023 

IF  (ABSIFM-FMH-. 00001)  13,12,12 

B 

25 

0024 

6 

IF  IA-l.OI  8,7,7 

B 

26 

0025 

. 7'  " 

FM-1.0 

B 

27 

0026 

B= ( 1 . +C2 1 

B 

28 

0027 

GO  TO  13 

B 

29 

0028 

8 

CONTINUE 

B 

30 

0029 

IF  (FM-l.OI  9,10,10 

B 

31 

0030 

. 2.... 

FMl-FM 

B 

32 

0031 

GO  TO  11 

B 

33 

10 

FM1*FM+0.1 

B 

34 

u 

CONTINUE 

B 

35 

0034 

GO  TO  5 

B 

36 

0035 

12 

FH1*FM 

B 

37 

0036 

GO  Til  5 

B 

38 

0037 

13 

POPI I I=B**1G/I i.-GJ  ) 

B 

39 

0038 

FMM( I I-FM 

B 

40 

0039 

14 

CONTINUE 

8 

41 

0040 

WRITE  16,151 

B 

42 

0041 

WRITE  17,161  INTG 

B 

43 

0042 

WRITE  17,171  (XRSI I 1 , POPI I 1 • I ,1=1, INTGI 

8 

44 

0043 

WRITE  (6,18)  IF  MM  I  1 1,  POPl'li  ,CST  Ak  ill  ,XKS<  II  .RMRSII) ,  I ,  I-l,  INTGI 

B 

45 

0044 

RFTURN 

8 

46 

C 

8 

47 

0045 

15 

FORMATI 1H1.2X, 'THE  DISPLACEMENT  THICKNESS  FROM  THIS  ITERATION  AND 

B 

48 

1THE  PRESSURE  DISTRIBUTION  FOR  THE  NEXT  ITERATION  FOLLOWS'//, 

8 

49 

28X,'M'tllXl'P/P0' ,11X,’0S'}13X, »X» ,12X, » RW» ,8X, • I  NTS'/) 

8 

50 

0046 

!<• 

FORMAT! I 31 

B 

51 

0047 

17 

FORMAT! 1P2E12,6,53X,I 31 

a 

52 

0046 

18 

FORMAT! 1P5E14,5,|6) 

TT 

53 

0049 

END 

B 

54* 

51 
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FORTRAN  IV  G 

LEVEL 

20  BEGIN  OATE  -  722B6  22/52/18 

OOOl 

SUBROUTINE  BEGIN 

C 

1 

0002 

COMMON  /GEN/  P£I,AHI,AHE,DP0X,PREF(2) .PR (2) ,P (21 ,OEN, AMU, XU, XD, XP, 

C 

2 

1XL,DX, INTO, CSALFA, ALPHA, XR.REORS, CAM, ZETA, PPO, TWTO, YSTART, USUP, 101 

C 

3 

2MEN, lHEAT,Z*TU*XSTEP/I/h,NPl«NP2»NP3»NEQ»NPH,KEX*KIN»KASE,KRA0/B/B 

C 

4 

• 

"3E  TA  iSaMAI  2 1  *  TALI  I  *  TAUC  *  A  jl  ( 2  )  *  AJE  (2 1  •  I  NO!  1 2ft  I NDE( 2  l/V/U  (20(F)  >  FI TTi 

C 

5 

4001, R C 200 )* RHOC 200 >* OMC 200 ).Y(200> 

C 

6 

C 

PROBLEM  SPECIFICATIUN 

C 

7  " 

0003 

REAO  (5,281  KRAD,IDIMEN,NEO,KEX,KIH, IHEAT.N 

C 

8 

000* 

READ'  (5,29)  REORS,  ZETA,  P  R  ( 1 »  ,  GAM,  ALPHA,  XR,  XL,  Us"UP,YST"ART,TWT6,X'STE" 

C 

9 

t 

IP 

C 

10 

c 

I DIMEN-0  FOR  PLANAR  FLOW  ANO  I0IHEN«1  FOR  AXISYMMETRIC  FLOW 

C 

11 

c 

INITIAL  EDGE  OF  BOUNDARY  LAYER  IS  YSTAR7 

C 

12 

0005 

PREF(1J-PR(1) 

C 

13 

c 

APPROXIMATE  CALCULATION  OF  UFDGE  FROM  ONE  DIMENSION  FLOW  RELATIONS 

C 

14 

0006 

UUMXX-S0P.T(2,/(GAM-i, I l*( (GAM*1. i/2. )**( IGAM+1. ) / (2.«(GAM-1. 1)1 

C 

15  ’ 

0007 

UUMXX-(1./(1.4XR) )**( 1+I0IMEN)/UUMXX 

C 

16 

0008 

UEDGE-UUMXX*SQRT(2. 1 

C 

17 

0007 

KASE-2 

C 

18 

0010 

I  F  I K  I N. E 0. "1 . OR  . KE  X, E 0 .  1 )  KASE*  1 

C 

19 

0011 

xu-o. 

C 

20 

0012 

NPH-NEQ-1 

c  " 

21 

0013 

NP1-N+1 

c 

22 

0014 

NP2-N+2 

c 

23 

0015 

NP3-NO 

c 

24 

c*  ■ 

INITIAL  VELOCITY  PROFILE 

c 

25 

0016 

Y! 1 1*0,0 

c 

26 

0017 

U(l)-0.0 

~tm 

27 

0018 

XNP2*NP2 

c 

26 

001* 

DEL Y* YST  ART/XNP2 

t r 

29 

0020 

00  1  1-2, NP3 

c 

30 

0021 

V(II*Y(I-1»»0ELV 

c 

'31 

0022 

I 

ETA-Y(I)/(0ELY*XNP2) 

c 

32 

0023 

U(l)-(2.*iTA-EfA**'2)*UEDGf 

c 

33" 

0024 

1 

CONTINUE 

c 

34 

c 

CALCULATION  OF  SLIP  VELOCITIES  AND  DISTANCES 

TT 

35 

0025 

BETA-1.0 

c 

36 

0026 

GO  TO  (2,3,4),  KIN 

c 

37 

0027 

2 

U(2)-U(3)/(l.t2.*BETA) 

c 

38 

0028 

Y 12 l»Y I  3) -BETA/I  2. ♦BETA) 

c 

39 

0029 

GO  TO  6 

c 

40 

0030 

3 

UU-U(1)*U(1) 

~c~ 

4l~ 

0031 

UI3-UI 1)*U(3) 

c 

42 

0032 

U33-U( 3) *U(  3) 

c 

43 

0033 

SQ»84.*U11-12.*U13*9. *U33 

c 

44 

.  0034 

U(2l-<  16.  •Uil-4.*Ul>U33)7(2.*fUl l)«U»f  1  ♦  SQRtTS'fr)  f 

■f 

V}  " 

0035 

Y(2)-YI3)*(U(2)4U(3I“2.*U(1))*.5/(U( 2)*U(3)4U( 1 1 ) 

C 

46 

0036 

GO  TO  5 

c 

47 

0037 

4 

IF  (KRAD.NE.O)  GO  TO  5 

c 

48 

0038 

U(2)-(4.*U(1)-U( 3)1/3. 

c 

49 

0039 

Y(2)-0. 

c 

50 

0040 

GO  TO  6 

“c" 

51" 

0041 

5 

U(2)-U( 1 ) 

c 

52 

0042 

Y(2')-Y(3T7~3. 

~c~ 

'53 

0043 

6 

GO  TO  (7,8,9),  KEX 

c 

54 

0044 

7 

U(NP2 1  -LMNPi )  /  ( i.*2.*BETA( 

c 

'55'” 

0045 

Y(NP2)«Y(Nj»3)-(Y(NP3)-Y(NPl))*BETA/(2.4BETA) 

c 

56 

52 
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00<-6 

GO  TO  10 

C 

57 

0047 

8 

U 1 1“  U  ( NP  1 1  *U  (  NP  1 1 

C 

58 

0048 

U13»UINPll*UINP3» 

C 

59 

0049 

U33=U(NP3)*U<NP3» 

C 

6b 

0050 

• 

S0=H4.*U33-12.*U13+9. *UL1 

C 

61 

0051 

U(NP21*( 16.*U33-4.*U13+U11»/12.*(U(NP1»*U(NP3> (♦SORT  (SQM 

C 

62 

0052 

Y ( NP  2 1 =Y( NP J 1 - ( Y( NP3) -Y(NP 1 N *( U( NP2» ♦U  ( NPlV -2i«>U ( NP3I  )  *.5/C0lNP2V  ' 

C 

63'""" 

1+U( NP  1 1  ♦  U (  NP  3 1  1 

C 

64 

0053 

GO  Tl)  10 

C 

65 

■  0054 

9  f 

U(NP2l  =  (4.*U(NP3I-U(NPlM/3* 

C 

66 

0055 

Y(NP2I«Y(NP3» 

TT 

67 

0056 

10 

CONTINUE 

C 

68 

0057 

IF  (NEO.EQ.ll  CO  TO  20 

C 

69 

0058 

00  19  J-ltNPH 

C 

70 

C 

INITIAL  PROFILES  OF  OTHER  DEPENDENT  VAR'lA'BlTS 

C 

71 

0059 

DU  11  1=1 »NP3 

c 

72 

0060 

ETA»YUI/(DELV*XNP2» 

c 

73 

0061 

F(J, II«THT0*(2.*ETA-ETA**2>*( l.-TWTOl 

c 

74 

0062 

11 

CONTINUE 

c 

75 

C 

CALCULATION  OF  CORR ESP0N3 1 NG  SLIP  VALUES 

c 

76 

0063 

GAMA ( J 1 = 1 .0 

c 

77  ' 

0064 

GO  TO  (12.13,141,  KIN 

c 

78 

0045 

12 

F  ( J ,  2 1  «F  (  J  ,  l" )  ♦  (  F  ( j  ,3 1  -F  ( J  ,  1  It  *(  1.  *BE  T  A-GAHAUH/ (  1 .  *BET  A*GAMA( J )  ) 

c 

79 

0066 

CO  TO  15 

c 

80 

0067 

13° 

G=(UI  2I»U(3I-8.*U(  111/(5.  *(U(  2  >*U(3»  U8.*UI1M 

c 

81 

0068 

GF  =  (  l.-PREF ( J 1 1 / ( l.*PREF(  J)  I 

c 

82 

006° 

GF«(G+GF)/| l.+G*GF I 

c 

83 

0070 

F(J,2I=F( J,3I*GF*(l.-GF»*f(J,l» 

c 

84 

0071 

CO  TO  15 

c 

85 

0072 

. IV... 

F(J,2I-F( J, 11 

c 

86 

0073 

IF  (KRAD.EQ.Ol  F ( J , 21 « ( 4. *F ( J , l> -F (J , 3 II / 3. 

c 

87 

0074 

15 

GO  TO  (16,17,181,  KEX 

c 

88 

0075 

16 

F ( J , NP  2 1 - F ( J , NP  3 1 ♦ ( F ( J  ,  NP 1 i -F ( j ,NP3»  I  •  ( 1 . ♦BETA-GAMA( j 1 J/( 1 . ♦BET AVG' 

c 

"89 

1AMAI J  >  1 

c 

90 

0076 

GO  TO  19 

c 

91 

0077 

17 

G=(IM  NP2)  *U(NP  1 1-0*  *U(NP3  11/ ( 5.  *(U( NP2I  ♦  U(NP1II*8»*U(NP3H 

c 

92 

0078 

GF=( l.-PREF( Jl l/( 1 . +PREF ( J 1 1 

c 

9  3 

0079 

GF-IG+GFI/I 1,»G*GFI 

c 

94 

0060 

F(J,NP2I=F(J,NP1 I*GF+ ( 1,-GF l*F( J  »NP3 1 

c 

95 

0081 

GO  TO  19 

c 

96 

008  2 

18 

F(J,NP2I»(4.*F( J,NP3I-FU,NP1I 1/3. 

c 

97 

00?  3 

19 

CONTINUE 

c 

9B 

0084 

20 

CONTINUE 

c 

99 

009  5 

CALL  DENSTY 

c 

1Q0 

C 

CALCULATION  OF  RAOI I 

c 

101 

0086 

CALL  RAO  (XU.RIll iCSALFA) 

c 

102 

0057 

IF  (CSALFA. EQ.O.. OR. KRAD.EQ.Ol  GO  TO  22 

c 

103 

0083 

PO  21  1=2, NP3 

c 

104 

0059 

21 

R(I»=R(1I-Y(  I>*CSALFA 

c 

105 

C 

CHANGE  HADE  IN  STATEMENT  NUMBER  28  FOR  INTERNAL  FLOW 

c 

106 

0090 

GO  TO  24 

c 

107 

0091 

22 

00  23  1=2, NP3 

c 

108 

0092 

23 

R ( I 1 =R (11 

c 

109 

0093 

24 

CONTINUE 

c 

110 

c 

CALCULATION  OF  OMEGA  VALUES 

c 

111 

0094 

0M(1I«0, 

c 

112 

53 
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0095 

I3MI2I-0. 

C  113 

0096 

DO  25  1*3. NP2 

C  114 

0097 

25 

ciN(ii>om  i- i i ♦. 5*1  rhoi  i  I»sm I i VRHorf- ii*ori-f ) WT i-i  n*TY  m-V' ■* c  ns*'" 

11 1-1)1 

C  116 

0098 

PEI-0M1NP2I 

C  117 

0099 

DO  26  1*3, NP1 

C  118 

0100 

26 

0Hin*aNm/PEi 

C  119 

0101 

0MINP2I-1. 

C  120 

0102 

0HINP3I-1. 

C  121 

0103 

IF  lflEQ.EQ.il  RETURN 

C  122 

**  o  104 

00  27  J* 1 . NPH 

C  123 

0105 

IF  (KEX.EQ.ll  INDECJI-1 

C  124 

0106 

IF  IKIN.EQ.il  INOUJI-1 

C  125 

0107 

27 

CONTINUE 

C  126 

0108 

RETURN 

C  127 

C 

C  128 

0109 

28 

FORMAT | 61 1 . I 31 

C“l29 

0110 

29 

FORMAT  (8E10.0) 

C  130 

0111 

END 

*C  131-*" 
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COSFF 


OAT E  *  72286 


22/52/18 


0001  SUBRUUTINF  CDEFF  0  I 

0002  COMMON  /GEN/  PE  I , AMI , AME ,DPDX, PRFF 121 1 PR (21 , P 1 21 .DEN, AMI'XUt XD, XP*  0  2 

. *IXL ,  OX,  IN  TG,CS  ALFA ,  ALPHA  ,XP *  RE OKS ,  GAM,£ETA,PPO,  TWTO, YSTART  «USUP, 1DI  D  3 

_ ZMEN,  lHEAT,Z,TO,XSTEP/l/N,NPl,NP2,NP3,NE(i,NPH,KEX,KIN,KASE,i(RAC/6/B  0 _ 4 

3ETAVGAMA(  2)  i  TAUI  ,T  AUE  •  AJI  ( 2 1  TaJEIZ I  ,  I ND I  (  2 1  ,~l  NDE  (  2 1  / V/U I 200 1  ,F(2.2"'  6  5 

*  400 I ,R I  2001 <RHO( 2001  ,0M( 200) i Y(200) /C/ SC 1 2001 » AUI200) t  BUI 200 1 , CUI 20  0  6 

. 50  )  t  A  (  2 r  200 1 1  B  (  2  >  200 1  ,  C  (2. 200 ) . . . D  7* 

0003  COMMON  /L/  AK.ALMG  _  06 

0004  .  DIMENSION  GH200I.  G2<200 1 1 '  G3<  206V."  6(  2 1 2'00i  •  $11266),  S2<200»,  $  D  9 

D  10 


• 

C 

CALCULATION  OF  SMALL  C  *S 

0 

11 

0005 

00  1  1-2, NP1 

D 

12 

0006 

RA«.5*IRI  I+ll+RI I  II 

0 

13 

C007 

KH«.5*(RH0I  I+1)+«H0(I)I 

0 

14 

0008 

UM=.5*IUII+1)+UII>) 

0 

15 

0009 

CALL  VEFF  (1,1+1, THUI 

D 

16 

0010 

1 

SCI  I )  =R A*kA*RH*UM*E  MU/ 1  PE  I *PE  1 1 

D 

17 

C 

THF  CONVECTION  TERM 

D 

18 

0011 

SA-RI 11*  AM  I/PE  I 

0 

19 

0012 

SB-1  R 1 NP  3 1* AME-R I 1 i *AMI 1 /PEI 

0 

20 

0013 

OX- xo-xo 

0 

21 

0014 

on  4  1-3, NPl 

D 

22 

0015 

GhD-UNd  +  ll-OMI  I-1I 

D 

23 

0016 

P2-. 25/OX 

D 

24 

0017 

P3-P2/0M0 

D 

25 

0019 

Pl*(OM(  I+ll-0MIin*P3 

D 

26 

0019 

P3*( DM! I l-OM( I-l 1 l*P3 

D 

27 

0020 

P2-3.*P2 

0 

28 

0021 

o-sa/ohd 

0 

29 

0322 

R2— SB*. 25 

D 

30 

0023 

R3-P2/UM0 

0 

31 

0024 

Rl—IUHI  1  +  11+3. *OM(II  1*R3  . 

0 

32 

0025 

P.3=(0M(  I-l  1+3.*0M(  l  ||*R3 

D 

33 

0026 

GUI  I-P1+Q+R1 

D 

34 

’  0027 

G2II 1-P2+R2 

D 

35 

0028 

G3(I  1-P3-0+K3 

0 

3b 

0029 

CUII1—  P1*UI  l+ll-P7*U(l  )-P3*UI  I -11 

D 

37 

C 

THE  DIFFUSION  TERM 

D 

38 

0030 

AUC I 1-2. /OMO 

D 

39 

0031 

BUI  I l=SC( I-1)*AU( I l/IOMIl »— OM <1  —  11 1 

D 

40 

0032 

Aum-scc  n*Auii  )/(cm<  i+u-uMii)  i 

D 

41 

0033 

IF  I NEQ.EQ. 1 1  GO  TU  3 

0 

42 

0034 

GO  2  J-i.NPH 

D 

43 

0035 

CIJ, I 1— Pl+FI J,I+1»-P2*FC J,1 )-P3*F(J, I-l) 

D 

44 

0036 

CALL  SOURCE  <J,1,CS,D| J,I 11 

D 

45 

0037 

Cl  J, I i*-C( J, 1 1 +CS-F 1  J,  I  l*DI J, I) 

D 

46 

0039 

A ( J, I  l-AUI 1 1/PREFI  j) 

D 

47 

0039 

BIJ, I l-BUl I l/PREFI J) 

D 

48 

0040 

2 

CONTINUE 

0 

49 

C 

SOURCE  TERM  FOR  VELOCITY  EOUATION 

0 

50 

0041 

3 

Sll I»=OPDX»OX 

D 

51 

0042 

S2I!I»P2*S1I  Il/IRHOII  )*UI  I  11 

0 

52 

0043 

S3! I )-P3*Si( I)/«RHU«I-1 l*UI I-l) 1 

0 

53 

0044 

S1II)*P1*S1<1 1 /( RHOI I+11  +  UII+ll  I 

0 

54 

0045 

cut  i  i*-cu<  n-2.*isi(i  )  +  s2<  II  +  S3HI  i 

0 

55 

0046 

snn-sii  n/uii+i) 

0 

56 

55 
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0047 

< 

S2(  1  )«S2(  n/UCl) 

D 

57 

0048 

S3(ii-S3(  u/uu-n 

0 

5B 

0049 

4 

CONTINUE 

0 

C 

COEFFICIENTS  IN  THE  FINAL  FORM 

D 

60 

0050 

00  5  I-3.NP1 

0 

61 

0051 

RL-1./(G2(I  l«AU(ll*BU(l)-S2(IH 

D 

62 

0092 

Aum-(Aumtsim-Gimi*RL 

0 

63 

0053 

BU(1I-(BU( I >  *S3( I )-G3( III 4RL 

0 

64 

0054 

5  . 

CU(II«CU(II*RL 

0 

65 

0055 

IF  (NEU.EU.1I  GO  TO  7 

0 

66 

0054 

00  6  J-l.NPH 

0 

67 

0057 

00  6  1-3, NP1 

D 

6B 

0358 

RL-1./(G2(1)*A( J.l l+B(J,I )-DCJ,I)» 

D 

69 

0059 

A 1 J , I t-IAl J, I  1 -Git  I >1*RL 

0 

70 

0040 

B(J,I>-(»< J.II-G3I ll»*RL 

0 

71 

0061 

6 

C  ( J 1  I  )  »C  (  J  >  I  1  *R  L 

0 

72 

0062 

7 

CALL  SLIP 

0 

73 

0063 

RETURN 

D 

74 

0054 

END 

0 

75- 
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FORTRAN  IV  C  LEVEL  20 


REAOV 


DATE  -  72286 


22/52/18 


0001  SUBROUTINE  READY  E  1 

0002  COMMON  /GEN/  PEI .AMI , AME, OPDX .PREFC2 I ,PR(2» ,PC2) ,DEN, AMU.XU.XO, XP,  E  2 

. "lXL.DX, IN  TG , C  S AL  FA , ALPHA ,  XR  ,' RE  ORSVGAM  VIETAVPPOYTHT  0, YST  ART , USUP, 1 0 1 '  E"  3' 

2MENi I  HEAT, I ,  TU,  XSTEP/V/UI 2001  ,F 12,2001 ,R (2001  ,RHO(200 1 ,0H(200 > , Y( 2  E  A 
3001/1  /  N,N  P 1VNP27NP3  ,"NEU,NPH,  KEX  ,  K1 N,  KASE  ,  K  ft  AO/B  /  B  ETTk  ,"  GAM  A  C  21 ,'  TAUT",  6  5  ' 

'  ATAUE , A J I ( 21 , A JE ( 2 1 , IN01 (2 ) ,  I  NDE 1 21  E  6 

CALL  DENSTY  . E  '“  7 

CALL  RAD  I XU,R( 1 I , CSALFA)  E  8 

. C . Y  NEAR  "THE'  I  BOUNDARY* . ** .  f  9 

_ •  GO  TO  C 1 , 2, _3I_,_K|N _ _ E _ 10 

1  VC  21 » (  f.  *BE_f  A  I  *0MC3I*4./C  C 3. *RH0(2 I ♦RHOI 3 1  I *CUI2 1 *UI3> 1 i  I  lT 

GO  TO  A  E  12 

. 2“ . Y(  2 )» 12. *  0MC"3 1/ ( YB'.'iR H0( 2  ) rRH 0 (3 ri* ( U1 2 >  VuCY|74."*U C  l > V f . E'l3~ 

GO  TO  A  E  1A 

. 3 .  Y(2l*.5*OM(3l/‘{RHOIll*UU7i . E  15 

A  YC 3 1 ■ YC  2 ) ♦•25*0M( 3 1 *( 1. /(RH0I3I *U(3I >42* /(RH0C3) *U(3I +RH0C2I *U( 21  I  E  16 

- iT  r~ir 

C  Y  'S  FOR  INTERMEDIATE  GRID  POINTS  E  18 

"00l"2* . DO  5  l-A.NPl .  . . "E  "19 

0013  5  VCII-YC I-1I*.5*(UM(I l-OMC I-l|)*(l. /(RHO(II*U( II 1*1. /(RHOI I-1I*U( I-  E  20 

. Kim . . e  "2i" 

_ C  Y  NFAR  THE  E  BOUNDARY _ _  E  22 

0014  YC  NP2 1 -Yl I  NPl  lV". 25* COM C NP2  I  -GM fNplII  *( 1./ CRHDCNP1I  *UC NP1 1 1  ♦2'.'? ( RHOT  E  23~ 


0003 

0004 

.  0005_ 
0006 
0007 
0008 
0009 
0010" 
0011 


0015 

0016  b 


YC  NP2 1  -Yl I  NPll  V".  25*C0M  C  NP2  I  -OMCNplII  *U./ IRHDINP1I  *UC  NP1 1 1  ♦2V/IRHS1*- E 
1NP 1 l*U(NP 1 1 *KHO( NP2 l*U( NP2I I )  _  E 

"GO  TO  (6,7,8),  KEX  . E" 

YCNP3I-YCNP2IM 1.+BETAI*(0MCNP2I-0M(NP1J l*4./((RH0(NPi)*3.*RH0(NP2  E 

1M*(U(NP1I»U(NP2I  I  I  ‘  **" . * .  E 

GO  TO  9  E 


2A 

25 

26 

27 

28 
29* 
30 
31* 


0017 _  _  . _ _ _  __ 

0018  7  Y (NP3 1 « Y ( NP  2 )  Tl  2 . * ( OM ( NP2 } -OMlNP 1 1) / C (RHOINPl I *3. *RHO INP21T*TTJTnP2  E 

lltUCMPll»4.*U(NP3M  I  E 

0019 .  GO  TO  9  .  . . E 

0020  %  Y(NP3l*V(NP2l*«5*(UM(NP2l"0M(NPlll/ (RHOI NP3I *U(NP3 I )  E  32 

"  0021 . 9* . IF  (  C  SAL  FA.  EO.  bi'.dR.KKAOi  EQ.'0l"'G'O''tO"l  l" . E'"33 

_ C  XXX  IS  USED  TO  KILL  SHOT  IF  NECESSARY _ E  3A 

0022  00  10  1-2, NP3  E  35 

0023  XXX-RC 1I*R(1I-2.*YCII*PEI*CSALFA  E  36 

0024* . IF  ( XXX, L T. 6. 0 f ' X~D-2~V*XL . E  *"37 

0025  IF  (XXX.LT.O. 01  GO  TO  1A  E  38 


0026 


To"  YC  1 1  »2.*Y(  1 1  *  PEI  /  (R  (  UVSQATCRC  1 1  <1  V-iV  *'Y  ( I  l*PE  I  *CS  ALF  A I )" 

C _ CHANGED  SIGN  OF  2  IN  THE  DENOMINATOR  OF  ABOVE  FOR  INTERNAL 


FLOW 


0027 

0028 

0029" 

0030 

0031 


11 

12 

13 


GO  TO  13 
DO  12  1-2, NP3 
"yc  n.PEi«Ycn/Rcrf 
Y(2I=2.*Y(2)-Y(3I 
"YCNP2J-2.*YCNP2l-YCNPli" 

CALCULATION  OF  RADII  _ 

I-"2,NP3 


0032  DO  14 

0033  IF  (KRAD.EO.OI  R ( I  I «R ( 1 1 

003A . IF  (KRAD.NE.0Y  RCI  l-RUl-YI  lT*CSAL'fiT . 

C  CHANGED  SIGN  IN  EXPRESSION  ABOVE  FOR  INTERNAL  FLOW 

00  35" . 14  C  ONT I NUE . . 

0036 _ IF_  C  R  CNJP3J  .LE._0.  01  XD»2.*XL _ 

CYCNP3f.~Lf.0V0l  XD"-2  ."•XL 


0037 

0038 

0039 


IF 
RETURN 

"end . 
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AEDC-TR-73-52 


FORTRAN  IV  G  LEVEL  20  OENSTV  DATE  -  72296  22/52718 


0001  ‘  SUBROUTINE  OENSTV  .  "  '  f '  ‘  1  ' 

0002  COMMON  /GEN/  PE  I .  AMI ,  AME,  O'PDX ,  PREF « 2 » , PR  1 21 ,  P  <2 1  ,OEN,  AKU,XU,XO,  XP,  F  2 

. "iXL.DX,  INTG tCSALFAt ALPHA »XRt  RE  ORSYGAmVZETA*  PPL)  i  TWT  0*  Y  START >USUP»[ D1  F  ’3," 

_ 2MFN, lHEAr,Z,TO,XSTEP/V/U<2001 ,F I  2 ,200 ) . R 1200 ) ,RH0(200 > , QM1200 1 , VI2  F  4_ 

*00l/I/N,'NPl,'NP'2,NP3',NEa,'NPM,KEX,'KfN,KAS¥7KRAT  F  5 

0003  v TNP3»Fll,NP3l-U(NP3l«*2/2.  _  F  6 

"  0004 . >H0NP3*TNP'3**fl.'/IGAM-iVi'r . " . 7  ' 

0005  DO  1  I-l.NPj  F  9 

ooot .  t»f<i,  n-u< i . f""9"‘ 

_ ;0007 _ 1  R  HD  I  I  |«RH0NP3«TNP3/T _ F  10 

0008  RETURN  F  11 

0009  END  _  F  12- 


FGRTRAN  IV  G  LEVEL  *20~  ~  "ENTRN  DATE  »  72286  “  '""22/ 52/1*9 


0001  .  SUQKOUTlNE  ENTKN  .  G  “  1 


. 1XL.DX, 1NTG, CSALFA, ALPHA, XR.RECRS, GAM, ZETA,PPO,TwTU,YSTART,USUP. 101  G  3 

_  _ 2kfn.IhFAT.Z»TU,XSTFP/V/U(200)  .F  ( 2  .200  I , R 1200 >  ,RHO( 200 1 1 0HI20Q > .  Y( 2  G  *  

3:  0  )  /  I  /N,NP1,NP2  ,NP3,NE0,NPH,KEX,  kTn.KASE  ,KRA0  "  *  '  G  ~  '5 

0003  1  GO  Til  ( 2i 3<6 1  ■  KEX  G  6 

0004  . 2  . RETURN  G  7 

0005  3  CONTINUE  G  8 

. C . THE  FOLLOWING  AHE  IS  FOR'  LAMINAR  'Flow .  G  "  9 

_■ _ 0006  _ IF  <INTG.NE.11  GO  TO  5 _ G  10 _ 

0007  00*4  l» 1 i NP3  G  11 

0009  IF  DMU  1.GT.0.9I  N9-I  G  12 

000°  .  IF  ISMI I  I.GT.0.91  I-NP3 .  *  G  13 

0010  4  CONTINUE  G  14 

0011 .  DaM»;)Nmi-0MiN9-u  . .  G  15"'" 

_ 0012 _ 0C'M9*0MIN9  1-3.9 _ G _ 16 

0013  5  CONTINUE  .  G  17 

0014  U9*Ul N9I -0UM9/UUM* I U( N9 I -UI N9-1 1 1  G  18 

0015  . TERM-3*  <U(N7I-UIN9-1  II /U’1  G  19"" 

0016  K9sKIN9l-r/T19/00'**IK|N9)-R<N9-ll)  G  20 

0017 . RHG9«rtH(HN9l-0CM9/UG-1*I  HUOI N9 1 -RHOI N9-1 1 V .  G  21 

_ OOlfl _ VIS9=VISCIUN9I-DI;MV/OU'|:*IVISCOIN9I-VI  SC0IN9-111 _ G _ 22 _ 

"0019  CUUP-P  It:91*RIN9m'H0IN9i"*UlN9l*VISC0fN91  G  23 

0020  CUUP=CUUP*R9*R9*RHfJ9*U9*VlS9  G  24 

0021  . •" . CUUP=CUUP/t2.*PFI  1  .  G"  25 

0022  CUU^*RI'M-ll*RIN9-ll*RH0IN9-ll*UIN9-ll*VISC0IN9-l)  G  26 

0023  "  . CUUtf=CUU-l*R<>*R9»‘iHU9*US*VlS9  G  27 

0024 _ CUUM»CUU'‘./I2.*PEll  _  _ _ G_  28 _ 

"“'0025  G5=2.*CU'JP/I0UM*!)UM9I  _  G  29 

0026  06=2. *CUUM/I OQM*I 0. 9-GMI N9- 1111  G  30 

0027  .  TERMA  =  G5* I UIN91 -U9 1-G6A IU9-UI N9-1 1 1  G  31 

002R  >  TEP.MA=TERMA/TEKMB  G  32 

■'  0029  . U9WIG»US'JP*SURTI2.  I  «S0RTI  U-PPO«*Y I  GAM^lYl'/GAMl"! . G"'33‘" 

C  USUP  IS  READ- IN  IN  BCGI_N _ IT  SUPPRESSES  THE  B.L. _ G  34 

0030’  TFRMCS IU9W IG-US  l/0X*0PDX/I  RHU9*U9l  *G  35' 

0031  TERMC=TERNC*PEI /TERMS  G  36 

0032 .  AME-TERrtA-TFRMC  *  . G"'"37*'" 

0033  AMF»AME-0.  l*Rl 1»*AMI  G  39 

00  J4 .  AME=AM6/I0.9*«INP311 . * . G  ”  39  '" 

_ 0035 _ RETURN _ G  40 

0036  6  AME'o'O.  G  4'f 

0037  RETURN  G  42 

0038  . END  . G  '43'-'"' 
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AEDC -TR-73-52 


FUKTRAn’i V~0  LEVEL-  ‘20  “  “rac  DATE  ■  72286  22/52/18 


OOOl 

SUBROUTINE  ESC  ( Xt J i 1 ND. A JFS) 

H 

1 

0002 

COMMON  /GEN/  PEI ,AM1 , AME i DPDX , PREF ( 2 1 ,PA (21 . P 121 «OENt AMU, XU .XD, 

XP, 

H 

2 

1XL , OX, INTO,  C SALFA, ALPHA ,X«,Rt ORS iGAM,2ETA,PP0,T«T0,YSTARTf USUP» 

ID1 

H 

3 

• 

2MEN , 1  HEAT ,2 »  TOtXSTEP 

H 

4 

c 

TW  IS  PRESCRIBED"  IF  1HEAT  »  1  —  UDDT  iSTFftFsClTiB'EDTI^fr'l 

H 

5 

000  3 

IN0«1 

H 

e> 

0004 

AJFS-TWTO 

H 

"7  " 

0005 

IF  (IHEAT.EQ.l)  GO  TQ  1 

H 

8 

-  0006 

IND-2 

H 

9 

0007 

AJTS-O.O 

H 

10 

0009 

1 

CONTINUE 

H 

11 

0009 

RETURN 

H 

12 

0010 

END 

H 

13- 

FORTRAN 

IV  G  LEVEL 

20 

NASS 

DATE  ■  72286 

22/52/18 

0001 

* 

SUBROUTINE 

NASS 

(XU'XD'AMl 

. i 

i 

C 

APPLICABLE 

TO  AN 

IMPF.RMEABL  E-WALL 

SITUATION 

i 

2 

0002 

AH“0» 

i 

3| 

0003 

RETURN 

i 

4 

0004 

ENO 

i 

5- 
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AEDC-TR -73-52 


FUR TRAN  IV  G  LEVEL  ’*20 


OUTPUT  "  DATE  -  722*86  *  ~  22/52/ 18 


0001 

0002. 


0003 
0004 
0005 
’  0006' 
0007 
*  0008 
0009_ 
0010 
0011 


0012 

0013 

_ 00 14_ 

00 1 5~ 
0016 
0017  ’ 
0015 
001° 
_  0020_ 
0021 
0022 
”  0023' 
0024 
0025 

_ 0026_ 

0027 
0023 
0029 
0030 
"  0031  * 

_ 0032_ 

0033 

...0034.. 

0035 

0036 

_ 0037_ 

003B 

0039 


0040 

0041 

0042 


SUBROUTINE  OUTPUT  J  1 

COMMON  /GEN/  PE  I , AM l , AME , DPDX.PREF 12 ) , PRI 2  I , P (2 ) , DEN, AMU.XU, XO, XP ,  J  2 
'1XL  ,PX,  IN TG, CSALFA,  ALPHA  ,Xk,REORS,' GAM,  ZE TA  ,PPU  ,TWTO, YSTART.USUP,  IDI  J  3[ 

2MfN,IHEAT,Z,TU,XSTEP/V/UI200)  ,F  (2 .200 )  ,R (200) , RHOI 200 > , OMI200) , LYI  2_J _ 4_ 

3001  /C/SCI200)',‘AU(  2001  ,*BUI  2bb)','CUl200)',A{2.200 1  .8(2  ,200  » ,012,200  I/O  j  *5 
4/YP ( 200)  ,UM200» , RM200I.HR I  2001  ,XM<200)  ,  PI  TOM  200)  .TEMPI 200 t/E/OS  J  6 
5TAR1300) ,  XRSC  3001,' RWRSi'300)  VC0SAL(300  )/  l  /N,NPl,NP2,NP3,  NEU,  NPH. K’6X*’*J  7* 
6,K IN, KASE .KRAD/B/BETA ,GAMA12 I , TAU1 ,TAUE , AJI ( 2 ) ,  AJE(2),IN0I(2), INDE  J  8 


7(2)  J  9 

_ IF  IlNTG.-NE.l)  GO  TO  1 _ J  10 _ 

A  L » AL  P HA  *"l  8  0.  >3'.  1 4 159 26*5  J  ~11 

WRITE  (6,7)  (0M(1),I-1,NP3I  J  12 

WRITt  (6,61  "KRAD,  1 D IHEN  .NEO.'KEXVKVN.THEAf,  N .  J*  13 

1  CONTINUE  J  14 

. DPDXG-OPDX*GAH/(GAM-I.  ) . * . J*  IS*'” 

IF  (KRAO.EQ.O)  DSTARI INTG I -Y( NP3 ) -PE  1 / (R ( 1 )*RH0INP3> *U (NP3)  )  J  16 

IF  (KRAD.'EQ.O)  GU-’  tO*  2  J  17 

DSTARI INTG)«R(l)-SQRT(R(l)**2-2.*CSALFA*(R(l)*Y(NP3)-0.5*CSALFA*Y(  J  18 

. 1NP3) **2-PE I/CRHCI NP3) *U(NP3I II )  .  J*  19*'" 

IF  (CSALFA. NE.O.Ot  DSTARI INTGI-OSTAR ( INTGI/CS ALFA  j'  20 

2  . "CONTINUE . .  *  J  ’”21’ 

_ IF  (FLO ATI  INTG-H/6..NE. FLOAT!  1  INTG-D/5))  RETURN _ J_22 _ 

DU* 3  !■ I.NP3*"  J  23 

TEMPI ll=FIl,II-UIII**2/2«  4  24 

. XM(  I  )»U(  I  )/SURT  ( (GAM-  1.  I*TEMP|  f  )")* .  J*  25*" 

XSQ«XM(I)**2  J  26 

.  IF  I XM(  l  I.LE.  1.  01  PITOT! I ) « (1  .♦(GAM-1.  )Vx*S0/2.'  )"** IGAM/I GAM-1.  ) ) . J  "  27 

_ IF  IXMIII.LE.1.01  GO  TO  3 _ J  28 _ 

P  ITOT ( 1 1 ■  (TgamVI .")  ♦  X S 0/2." I  **TG  A / ( G A M- 1 7T7  j 7  29 

P I TOT (  I ) «P I  TOT ( I ) * ( IG AM* 1 . I / ( 2. *GAH*X  SQ-GAH* 1 . ) )** 1 1 . / 1  GAM- 1, M  J  30 

*3* . CONTINUE  . j  31**" 

SQ2-SQRTI2. I  3  32 

DU  4  I-1.NP3  .  J  33 

_ P1T0M  IMPITOTt  H/PIT0TINP3) _ J_34 _ 

YRI 1 1 »Y(  (  l••CSALFA/R  1 1*1  "  J  35* 

UR (I )  «U( I)/U(NP3)  .  J  36 

. RR(I  I -RHOI  I  I  /kHO(  NP3I  " . *"' .  . J  37  '** 

HPII l-U( II/S02  J  38 

V"*  CONTINUE  .  .  .  J  39  * 

_ WRITE  !  6, 8 1 _ _  _  _ _ J _ **0  _ 

WRITE  1 6 , 12)  IN  TG,  XU,  l  ,RE0RS  ,'GAM ,  PR  ( IT7Z  E  TA  ,7  WTO  ,  AL, XR, XL  **J  41* 

WRITE  (6,131  NP3, DSTARI INTG), PE  I, AME , DPOXGiOX, CSALFA, USUP ,PPO .TAUX  J  42 

. *1,AJIU)  . * . . J*  43  * 

WRITE  (6,9).,.  J  44’ 

00  5  1«=1,NP3,"2  * . * . J  45  "" 

_ WRITE  (6,111  (Y(1),F(1,1),0K(1) ,HH1I)  ,R( I  I ,RHO( I) , YRI 1 ) ,XM( 1 ) .TEMP  J  46 

1(  I )  ,RR(  l")  , PITOT*(  I ) )  J  47 

5  CONTINUE  4  48 

. WPITE  (6,10)  '( Y"(NP3)  ,F(*iiNP3l  ,"UR(NP3*)','h'R(NP3")","R(NP3)  *,RH0'(NP3)',VR'(N  *  J*  49*" 

1P3),XM(NP3I,TEMP(NP3) ,RR(NP3) .PIT0TINP3) I  J  50 

(6,9) .  4  51 

RETURN  J  52 

C  ’J*  53 

6  FORM AT (////////, 20X, 'THE  INPUT  FLAGf  ARE •  ,/// ,25X, *KRAD  -  ',11,  J  54 

. *i//,25X,  MDIMEN  *-"»Vfi  ,//,25X,  'NEO . «  *' VU*,/*/ ,25X;  »KEX  V*  ' ,  IT,*  *  j"*55"" 

2//.25X, 'KIN  -  • ,I1,//,25X,’IHEAT  -  * , II ,// ,25X, »N  »  ',131  J  56 
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AEDC-TR-73-52 


FORTRAN  IV  C  LEVEL  20  “  OUTPUT  DATE  =  72286  22/52/18 


0043 

. 7"'*' 

FORMAT! 24H1T HE  VALUES  UF  OMEGA  ARE/I 1P10E 11.41 1 

J 

57 

0044 

. . B ..... 

FURMATIlHlfAXt’INTG'tOXt'XU’tilXt’l’iOXr 'REORS' »7X, 'GAMMA* , 8X, *  PR* 

J 

58 

1.9X. 'ZETA'.SX, 'TW/TO* ,7X, ' ALPHA • , 8X ,  • XR‘ ilDX. *XL*/,5X, 'NP3',9X, 'OS' 

J 

59 

2TAP 8X, 'PEI' »9X» *AME • , 8X , 'OP/DX » , 8X . • DX * , 8X , 'COS ALF » , 7X, 'USUP* , 8X 

J 

60 

3,'P/P0',8X,'TAUI',7X,  'AJl til1/) 

J 

61 

0045 

9 

FORMAKtiX,  'Y'.qX.'M/HO'.RX,  'UAIE'.SX,  'U/UM*  ,9X,*R',10X,  »RHO*  »10X, 

J 

62 

l'Y/R*,9X, '  M  * ■ 1 1 X  r ' T ' • lOXi »N/ME » ,7X, 'PITOT V) 

J 

63 

0046 

10 

FORMAT  1 1P11E12*5|/I 

J 

64 

.0047 

11 

FORMAT I1P11E12.5) 

J 

65 

0048 

12 

FORMAT 1 I8.4Xt 1P10E 12*4) 

J 

66 

0040 

13 

FORMAT  1 1 8 . 4X,  IP  10E 1 2.  4,//  ) 

J 

67 

005D 

END 

J 

68- 

FORTRAN  IV 

G  LEVEL 

20  PRE  DATE  -  72286  22/52/18 

0001 

0002 

SUBROUTINE  PRE  (X.OPOXX) 

COMMON  /GEN/  PE(,AMI,AHE,OPOX.PREF(2I ,PR(2) .P(2> ,DEN, AHU.XU.XD.XP, 
UL.DX, INTO, CSALFA, ALPHA , XR, REORS, GAM,ZETA,PPO,THTb,YSTART,USUP,IDr' 
2MEN, (HEAT. Z, TO, XSTEP/V/Ut 2001 ,F (2i200) >R(200I tRHO(200) , 0M(200),Y(2 

K 

K 

"K 

K 

1 

1 

• 

300)/I/N,NPl,NP2,NP3.NEQ,NPH,KEX,KlN,KASE,KRAD 

R 

5 

0003 

DIMENSION  XX (300) ,  P0P(300) 

K 

6 

0004 

IF  (INTG.NE.O!  GO  TO  1 

X 

"7" 

0005 

READ  15.3)  LMAX  .  . 

K 

8 

0006 

RE  AO  15.4)  UXILl.POPIU.L-l.UMfX) 

K 

9 

'  0007 

WRITE  (6,5) 

X 

10 

0008 

WRITE  (6,6)  (XX(L) .POPILl .L.L-l.LMAX) 

K 

11 

0009 

....JL _ 

CONTINUE 

K 

12 

0010 

L»1 

K 

13 

0011 

2 

CONTINUE 

K 

14 

0012 

L«L*1 

X 

15 

0013 

IF  (XX(L).LT.X)  GO  TO  2 

K 

16 

D014 

OPDX* (POP ( L) -POP (L-l 1 ) / (XX( L)-XX(L-1 ) 1 

K 

17 

0015 

PPO*POP( L-l) ♦OPOX*( X-XX(L-l)) 

K 

18 

0016 

DPDXXaOPDX 

K 

19 

0017 

DPDXX»0P0X*IGAH-1. I/GAM 

K 

20 

0018 

RETURN 

X 

21 

C 

X 

22 

0019 

3 

F0RMATI13I 

X 

?3 

0020 

FQRMATI2E12.0) 

K 

24 

0021 

5 

F0RMAT(1H1,/,8X,'X',12X,'P/P0',/) 

X 

25 

0022 

6 

FORMAT! 1P2E14*  5 . 1 6) 

X 

26 

0023 

END 

X 

27- 
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AEDC-TR-73-62 


FORTRAN  IV  G  LEVEL 

20  RAO  OATE  ■  72266 

"22/52/18 

0001 . V 

SUBROUTINE  RAO  I Xi R 1 i CALPMA I 

1  — " 

'i 

0002 

COMMON  /GEN/  PEI, AMI ,AME,DPDX,PREF<2I ,PH(2I ,P(2) ,OEN, AMU, XU, *0, XP, 

L 

2 

1XL,DX, INTG,CSALFA,ALPHA,XR,RE0RS,GAN,2ETA,PPU,TmT0,YSTAKT,USUP, 101 

L 

3 

2MEN, I  HEAT, Z  ,T'J«  XSTEP/V/Lit 2001 ,F(2,200I «R 12001 ,RHO(200I ,  I 

DM1200 1 ,  VI 2 

L 

V 

m 

3001/I/N,NP1,NP2,NP3,NEQ,NPH,KEX,MN,XASS.RR40 

L 

*  5 

c 

APPLICABLE  TO  NOZZLES  WITH  CONSTANT  LONGITUDINAL  RADIUS 

OF 

L 

6 

c 

CURVATURE  OF  THE  C0NVFPG1NG  SECTION-  CONSTANT  MALL  HALF 

ANGLE  OF 

L 

7  •" 

c 

DIVERGING  SECTION- A NO  HALL  SLOPES  MATCHED  DOWNSTREAM  OF 

THE  THROAT 

L 

8 

0003 

IF  IINTG. NE. 01  GO  TO  1 

L 

9 

•0004 

P 12-3. 14159265/2. 

L 

10 

0005 

ALPHA>ALPHA*PI2/90. 

L 

ii 

0094 

COSALF-COSt ALPHA  I 

L 

12 

0097 

S1NALF-SINIALPHAI 

L 

13 

onus 

ZWIG>XR*( l.+SINALFI 

L 

14 

0000 

XWIGaXK*( PI  2* ALPHA  1 

L 

15  " 

0019 

RHIGal.*XR*C 1 .-COSALF I 

L 

16 

0011  1 

CONTINOC 

L 

L7 

0012 

IF  (X.GE.XHIGI  GO  TO  2 

L 

18 

0013 

Rl  1)*1.*XR*(1.-SINIX/XR|| 

L 

19  '  ' 

0014 

CSALFA-SINIX/XRI 

L 

20 

0915 

CALPHAaCSALFA 

L 

21 

0016 

H  1-RI1I 

L 

22 

0017 

ZZ-XA*U. -CUSIX/XRN 

L 

23 

0018 

GO  TO  3 

L 

24 

0019  2  • 

CONTINUE 

L 

25 

0029 

P(1I"<X-XWIGI*S1NALF*RMIG 

L 

26 

0021 

CSALFA«COSALF 

L 

27 

0022 

CALPHAaCSALFA 

L 

28 

0023 

Pl-RUI 

L 

29 

0024 

2Z-(X-XHIGI*CSALFA*ZHIG 

L 

30 

0025  3 

CONTINUE 

L 

31 

0026  » 

IFUOlMEN.EO.OIRlal. 

0027 

Z2-ZZ-XR 

L 

32 

9028 

2-ZZ 

L 

33 

0029 

RETURN 

' 

L 

34 

0030 

END 

L 
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AEDC-TR-73-52 


FORTRAN  IV  G  LEVEL  20 


SLIP  DATE  ■  72286  22/52/18 


0001 

0002. 


00*03 


.  0004 

0005* 

0006 

0007 

■0009 

0009 

_OD10 

oo  ii* 
0012 
0013 
0014 
0015 

_ 0016 

0017 

0018 

001« 

0020 

‘"0021 

_ 0022 

0023 

0024 

0025 

0026 
_  *>027 
0028 
0029 
0030 
0031 
0032 
_0033 
0034 
0035 
"'0036 
0037 
0038 
—0039 
0040 
0041 
0042 

“6643 

_ 0044 

0045 

0046 

0047 

0048 


,  SUBROUTINE  SLIP  Hi 

Common  /gen/  pf i .ami ,ame,opdx,prem2> ,pri2>,p(2> ,den,amu,xu,xd,xp,  h  2 

. 1XL, nX.INTG.CSALFA, ALPHA, XR.REOR'sVGAH'.ZETA, PRO, THTO, YSTARTiUSUP,  IDI  H*  ‘  3 

2HEN,  IHEATtZ  , TQ»XST6P/  I /N_, NP_1^NP2 * NP3  ,NEO_jNPH,KEX,XIN»KASE,RRAO/V/U _ H _  4 _ 

"  3I200I  ,FI'2‘,200T.RI2*00»  ,RHO'l200)  .dHI200r,  VC26'0r/8/BETA,6AHAl'2r,  TAUT.  H  "5 


4TAUE, AJl I  21 , AJE42 } , INDI 12 ) , (NOE  121  H  6 

. r‘  COMMON  / L/”  A K » A L MG /C / SC "(2 0  0 1, AU 1 200* lYBUI 1 20& I » £ 0 1 2  00 l*»  Ai 2 , 2 0 0 )*, B f 2*,  '  M  7 

1200 1 ,C( 2,2001  M  8 

C .  SLIP  COEFFICIE'NiS‘NEAR"THET'BOW6ARV‘FOR'VELi»£i'fV  EUUAtTON . N  "  9‘ 

l  CUI2I-0.  _ M  _10 _ 

CU(NP2)-0”i  M  11 

GO  TO  1 1, 2,31  f  KIN  _  .  _N  12. 

Y"  '  BUI 2 1-0.  .  H  13 

AU(2I-1./(1.«2.*BETAI  M  14 

. GO  TO  5  .  M  15 

2  _ SQ-84.-UI 1)*UI  11-12. «U(  1I«U(3I*9.»LH3I*UI3> _ H  16 

BUI  2 1  ■  8  •  *  I 2."*Ul17*U(  3  1 1  /(  2.  *U  1 1 1  *7.  *U  ( 3 1  ♦SQRTiSQ)  I  H  it 

AUI 2 1 ■  1  •  —  liU I  21  H  18 

60  TO  5  .  '  *  M  19 

3  BUI  2 1-0.  M  20 

. CALL  VEFF  (2,3,FMUI*  '*  .  H  21*' 

.  AK 1- 1 . /DX-DPDX/I RhOI  1 1 »U1U*UU)  > _ H  22 _ 

AK2— Um*AKl*DPUX/(RHO(Il*UI  II  »  H  23 

AJ-RH0(LI«J(LI*.25«(Y(2)*Y(3)I**2/EHU  H  24 

V  IF  IKRAO.EO.OI  GO  TO  4  ‘  "  H  25 

AUI2)»2./(2.*AJ*AK1I  _  M  2ft 

*■  CUI2I— .5*AJ*AR2*AUI2I  "  *  .  . H"27"' 

_ GO  TO  5 _  _ H  28 

4  CUI2»-i./(2.*3;*AJ*AKll  H  29 

AUI2I«CUI2I*I2.-AJ*AKII  .  H  30... 

. CUI 2I--CUI 21 *4.*AJ*AK2  H  31 

C  SLIP  COEFFICIENTS  NEAR  THE  E  BOUNDARY  FOR  VELOCITY  EQUATION  H  32 

*5 . GO  TO  (6,7,81,  KEX  .  . M*'  33'" 

6  _ AUINP2I-0. _ _ M _ 34 _ 

BUINP2I-1./I 1.*2.*8ETaI  H  35 

GO  TO  9  .  N  36. 

7  . SO-84.  *U(  HP  3 1 -UINP3I  - 1 2.*  U I N  P  3 1  *U(N>1  i**;  YC'l'NP'i  fni'l  NP II . N"' 37. 

AUINP2 l-6.*( 2.-UINP3I «U( NPl I I /( 2. *U( NP3 I ♦7.*U INP1 I ♦SORT (SQ) I  H  38 

. BUI  NP 21-1. -AUINP2I  . M  39" 

_ GO  TO  9 _ H _ 40 _ 

8  AUI  NP 2 1-0.  H  41 

_  CALL  VEFF  (NP1,NP2,EMU)  H  42 

. BK  l-l./DX-OPDX/ 1  RHOI  NP3 I •UINP3I  -U(NP*3*i  f  . H"  43 

BK2— U(NP3l*8Kl*OPOX/ (RHOI NP3 )*UINP3I I  H  44 

. 8  j-RHOI NP 3 ) *U( NP  J)  *. 2 5*1 2 .  *YC  NP3I-V I NP  IT - X 1 N P~2 ) T* *2 / EMU . H* "  <5 

_ CUINP2I-  l./(2.*3j*8J*8Kl). _ M  46 _ _ 

BUI NP  2 1 -C  U( NP2 1  * 1 2. -8 J*6K 1 1  M  47 

CUI NP2I--CUI NP2I*4«*BJ*BK2  _  H  48 

9  IF  (NEQ.EQ.  1 1  RETURN  *'*  . M-  49  ' 

C  SLIP  COEFFICIENTS  NEAR  THE  I  BOUNDARY  FOR  OTHER  EQUATIONS  M  50 

'DO  20  J- 1  ,NPH  . H"  51  "■’ 

_ C I J ,  2  L»Ot _ H_  52_ 

CU.NP2I-0.  M  53 

. .GO  TO  1 1.0,12,  IAIjl.KIN. _ _  H  54„ 

10  CALL  FBC  I XD, J, INDI (J 1,011  H  51 

}F  I INDJ  ( J  r,EQAlJ  GO  TO  1*  N  56 


63 


AEDC-TR-73-52 


FORTRAN  IV  G  LEVEL  20  SLIP  DATE  -  72286  22/52/18 


0049 

AJII JI-QI 

H 

57 

0050 

A(J,2)«1. 

H 

5B 

0051 

BIJi 21-0. 

M 

59 

0052 

C(J,2)«8.*( l.*2.*BETA)*PREF< J )*AJI  (  J ) /< AK*AK*BETA*I l.+BETA»*l 1.+6E 

N 

6b 

ITA)*( 3,«RH0( 2 )*RHQ( 3) t*U<3t  1 

H 

61 

0053 

GO  TO  15 

H 

62 

0054 

ii 

FIJ,  11-01 

N 

63 

0055 

AIJ»2)-( 1.4BETA-GAHAI J > ) / 1 1 .*8E T A+GAMAI J I  I 

H 

64 

0056 

81 J,2)«l.-A< J,2) 

H 

65 

0057 

GO  TO  15 

N 

66 

0058 

12 

AIJ,2>»<UI2)4U(3I-8.*U( in/1  5.*l  UI2I+UI3  )  )*8.*U(  111 

N 

67 

005J 

GF«I  l.-PREF I J ))/(l.*PREFIJ)  ) 

N 

68 

0050 

A(J,2)«(  Al  J,2)+GF)/I 1  .♦  A 1 J ,  2  I  *GF  I 

N 

69 

00b  1 

8(J.2)>1.-Af J,2) 

N 

70 

0062 

GO  TO  15 

N 

71'  ' 

OOr.8 

13 

Bl J.21-0. 

M 

72 

0064 

CALL  SOURCE  IJ,1,CS,0S) 

H 

73 

0065 

AK1-1. /OX-OS 

M 

74 

0065 

AK2*-AKL«F< J.ll-CS 

H 

75 

0067 

AJF-AJ4PRFF1J! 

H 

76 

0943 

IF  IKRAO.cO.O)  GO  TO  14 

M 

77 

006° 

A< J,2l-2./(2.*AJF*AKLI 

N 

78 

0970 

C ( J. 2) 5*AJF*AX2*A| Jt 2) 

N 

79 

0071 

GO  TO  15 

M 

80 

0072 

14* 

C<J,2)>1./(2.»3.«AJF*AK1) 

M 

81 

0073 

AIJ.2I-CI J,2)*I2.-AJF*AK1) 

H 

82 

0074 

CIJ.2I— CIJ,2I*4.*AJF*AK2 

N 

83 

C 

SLIP  COEFFICIENTS  NEAR  THE  E  BOUNDARY  FOR  OTHER  EQUATIONS 

M. 

84 

0075 

15 

GO  TO  (15,13.19),  KFX 

M 

85 

0076 

16 

CALL  F8C  (XO,J,lNOE(Ji,OEI 

N 

86 

0077 

IF  I INDSI JI.EQ.ll  GO  TO  17 

N 

87 

0073 

s 

AJEI Jl-OE 

N 

86 

0979 

B(J,NP2)*1. 

N 

89' 

0080 

AI J.NP2I-0. 

N 

90 

0081 

CIJ,NP2l«-8.*t 1.42.*BETAI*PREF( JI*AJEI J)  /  C  AK*AK*BETA*I  1  .♦BETA  I*  ITT" 

N 

91 

1*RETAI*IRH0INP1)»3^*RH0(NP2) )*U(NP1) ) 

M 

92, 

0082 

GO  TO  20 

H 

93  ' 

0063 

17 

FI J.NP3I-0E 

H 

9h 

0034 

B(J,NP2)«C 1.4BETA-GANAI Jl ) / 1 1 .♦ BETA»GAM4(  J )  I 

M 

95 

0035 

Al J.NP2I-1.-BIJ.NP2I 

H 

96 

0036 

GO  TO  20 

H 

97 

0037 

18 

B(J.NP2!~<U(NP2)*U(NP1I-8.*UINP3M/(5.*(U(NP2I*U(NP1>  K8.*UtNP3l 1 

M 

98 

0083 

GF«I  l.-PRFF I J )  I  /(  l.»PREF(JI) 

H 

99 

0039 

BI J,NP?I=I Bl J,NP2)+GF)/(1.*B( J,NP2)*GF) 

N 

1Q0 

0090 

A(J,NR2)*1.-BIJ,NP2I 

H 

101 

00°1 

GO  TO  20 

M 

102 

0092 

19 

Al J,NP2)»0. 

H 

103 

0093 

CALL  SOURCE  ( J , NP3 ,CS |DS) 

N 

104 

0094 

8K1-1./9X-0S 

H 

105 

0095 

BK2»-BX1*F 1 J , NP3I-CS  ^ 

H 

10b 

0096 

BJF-8J4PREFIJI 

H 

107 

0097 

C  C  J , N P 2 1 ■ l./(2.«3.«BJF*HKl) 

M 

108 

0098 

0(J,NP2)-CU,NP2>*I2.-BJF*BIU) 

M 

109 

0099 

CIJ.NP2>~-C(J.NP2>*4.*BJF*8K2 

M 

110 

0100 

20 

CONTINUE 

H 

111 

RETURN 

H 

112 
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FORTRAN 

IV 

G 

LEVEL 

20  SLIP  DATE  -  72266  22/52/18 

0102 

END 

M 

iTS- 

FORTRAN 

IV 

G 

LEVEL 

20  SOLVE  OATE  -  722B6  22/52/18 

0001 

SUBROUTINE  SOLVE  { A,0  ,C fF ,NP3 1 

N 

i 

C 

THIS  SOLVES  EQUATIONS  OF  THE  FORM 

N 

2 

C 

Fin  ■  Ain«Fim>  ♦  bii)*fii-i>  ♦cm 

N 

,3 

C 

FOR  1-2.NP2 

N 

4 

0002 

DIMENSION  AINP3I  «  6INP3),  CINP3),  FINP3)  " 

N 

s 

0003 

NP2-NP3-1 

N 

6 

0004 

BI2I>BI2»*F(1)4CI2> 

N 

7 

000  S 

DU  1  1-3, NP2 

N 

8 

0006 

T-l./l l.-B(IJ*AII-lll 

N 

9 

0007 

All) -All) *T 

H 

10 

0008 

1 

BID-101  I )  *BI  1-1)4C(II)*T 

N 

U 

0009 

DO  2  1-2, NP2 

N 

12 

0010 

J-NP 2-142 

N 

13 

0011 

2 

F( J)-AIJ)*FI J41J48IJI 

N 

14 

0012 

RETURN 

N 

"IS 

0013 

END 

H 

16- 

FukThAN 

IV 

G 

LEVtL 

20  SOURCE  DATE  ■  73U4b  11/33/59 

0001 

SUBROUTINE  SOURCE  1 J, IfCSiDS) 

0002 _ 

.  COMMON  /GEN/  PEI,AMI*A«EiUP0XiPHEFl2» ,PR  ( 2  ).,P  (2) »DEN, A*U, XU, XU. XP, 

.A. 

...  3  ... 

1  XL, UX, INTG.CSALFA, ALPHA, AR, REORS, GAM, 2E  T  A,PPO, TwTOi YS I  ART, USUP, IOI 

A 

4 

_A_ 

* 

Jt  f  A  ,  RaH  A  (2)  ,  TAUI ,  TAUf_,  AJI  (2)  ,AJE(2>  fINDl  (2)  ,  INI)E  C2J/V/U  1200)  ,F  12,2 

A 

6 

„4 u 0 ) , HI20u),KHO(2UO),OM (200) ,Y  1200) /C/SC (200) t AO (200) ,UUt200) ,CUC20  _ 

A. 

..  7  .  .. 

50), A(2,2>»0) ,H (2,200) , C (2, 200) /O/YR (2JU) , UR (200) , HR (200) , HR (200) ,XM 

A 

8 

6(20  0)  tP I  TO I  (2i)U).,.TEMP  1200)  /E/DSTAR130D)  •  XliSX3QiU_«  KHRS  13.08.1 1.C0SAL.I3 

A 

.  H 

.700) 

A 

10 

onuJ 

0004 

CS-CS-SC  I  T-l ) • IU ( I ) *U ( I ) -U ( I-l ) *U ( 1-1) ) / I0M(I ) -OM 1 1- 1 ) ) 

000b 

CSa  11  .-1  r/PM EF  1  J)  1  *CS/  (Ort(  1*1  )  -OM  1 1-1  )  ) 

OOOe 

US-0. 

0007 

turn 

000a 

END 

FORTRAN 

IV 

G 

LEVEL 

20  VEFF  DATE  ■  72286  22/52/18 

0001 

SUBROUTINE  VEFF  (1,1 PI, EMU) 

P 

l 

.  -Q902..... 

— 

COMMON  /GEN/  PEI,AMI,AME,0PDX,PREF(2) , PR (2) ,P(2) , DEN, AMU, XU , XD, XP, 
IXL ,DX , INTO.C  SALFA , ALPHA ,XR ,REORS,GAM, ZtfA,PPO,TWTO,YST ART, USUP, IOI 

P 

-p- 

2 

— 

2MEN ,IHEAT,Z,TO,XSTrP/V/UI200l TF (2,200  I ,R (200) ,RHO(200) ,0M( 200),  VI  2 

P 

■4 

300)/l/N,NPl,NP2,NP3,NEQ,NPH,KEX,KIN,KASE,KRAO 

P 

5 

0003 

T—  F I  1, 1 l-UI I )**2/2, 

P 

6 

0004 

TT-FI 1 , IP 1)-U( IPl )**2/2« 

P 

7 

0005 

T-(T4TT>/2. 

P 

8 

0006 

EMU-T**ZETA/( REORS/SORTI2. ) ) 

P 

9 

0007 

RETURN 

P 

10 

0008 

END 

11- 
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FORTRAN  IV  C  LEVEL  20  VISCO  DATE  -  72286  22/82/18 


0001  . FUNCTION  VISCO  (I)  0  1 

0002  COMMON  /GEN/  PE  1 , AMI , AM£ , DPDX ,PREF 12  I ,PR (21 ,P (21 ,DEN, AMU, XU, XD, XP,  U  2 

. lXL.UX.lNTG.CSALFA, ALPHA, Xft , RF CRS, GAN, ZETA.PPO , TWTO.YSTART .USUP, ID1  0  ,3 

2KEII,  1  HFAT » 2  » TO,  X$T£  P/V/U(  200)  .F  1 2  ,200) .» (200 )  ,RHO( 200 ) ,  0MI20Q ) ,  YI 2  Q _ 4_ 

300l/(/N,NPl,NP2,NPJ  iNeU tNPHiKEX  t*KIN. K.ASE  t KRAO  U  5 

0003  T-F(l,I)-U< ll**2/2.  _  0  6 

0004  . V I  SCO«T**ZETA/(  REOR S/SQRTI 2.  ) )  0  7 

0005  RETURN  .  .  0  B 

0006' . END .  *  0  9- 


FOPTRAN  IV  G  LEVEL  20  WALL  GATE  -  72286  22/52/18 


0001 

...  0002  ... 

SUBROUTINE  WALL 

COMMON  /GFN/  PE  It AMl,AMf,0PDX,PREF(2l,PR(2>,P|2l,DEN,AMU,XU,XD,XP, 
1XL. CX,  INTG, CSALFA, ALPHA. XKt'RECkS. GAM, 2ETAtPPU,TWTO.VSTART.USuPtl6l' 
2WEN,IHF.\T,2,T0,XSTEP/V/U(200I  .F ( 2,200)  , R (200 ) .RHOI200 ) ,  OM  ( 2  00 1 ,  Y(  2 
3001/ 1  /  N, DPI  , NP2  .IJP3, NEO. NPH.ICEX.KIN.XASE  ,  KRAU/B/8E~TA,  GAMA(2  )  iTtAUl  , 
4TAUE.AJI  I2I.AJEI2),  INGH2I.IN0EI2) 

R  1 

R  2 

■  R-— ,3 

_* _ ,9 _ 

R  5 

R  6 

c  ■ 

C  ALCULAT  ID  i  OF  OFT  A  FOR  THE  I*  BOUNDARY 

R 

7 

0003 

1 

YI-.6«(Y(2t«Y(3H 

R 

B 

0006 

UI*.S«(0(2I»U(3II 

R 

.  0005 

• 

R  H». 25*1  3.*PH0( 2I«RH0I3)I 

R 

10 

0006 

RE-RH*Ul*Yl/VlSCO(ll 

R 

11 

0007 

FP=UPDX*YI/(KH*UI*UI) 

R 

12 

0008 

AMaAMI  /  ( RH*IJI ) 

R 

13 

c 

FOR  LAMINAR  FLOW  AND  AM-0  (NEED  DIFFERENT  EXPRESSION  IF  P-OI 

R 

14 

0009 

S-1./KE-FP/2. 

R 

.  i  5— . 

0010 

BETA«RE*(S*FP*AMI 

R 

16 

0011 

TAUl»S*RH*Ut*UI 

R 

17 

0012 

IF  (NEQ.Esl.l)  RETURN 

R 

18 

c 

CALCULATION  OF  GAMA  >S  FOR  THE  I  BOUNDARY 

R 

19 

0013 

00  2  J-l.NPH 

R 

20 

c 

FOR  LAMINAR  FLOW  ANO  LIMITING  ZERO  AM 

R 

21 

0014  . 

SF-l./(PA< J)*RF) 

R 

22 

0015 

GAMAI JI«RE*PR( JI*(SF*AM) 

R 

23 

0016 

IF  (INOHJI.EO.il  AJ1  ( JI*SF*RH*UI*(2»*F( J|i)-F( J»2)-FIJ,3II*.  5 

R 

24 

0017 

2a 

CONTINUE 

R 

25 

0018 

RETURN 

R 

26 

0019 

. 

END 

R 

27- 
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APPENDIX  IV 

INITIAL  PRESSURE  DISTRIBUTION  COMPUTER  CODE 


1:  C 
21 

3!  C 

4:  c 

5:  C 

5:  c 
c 

s:  c 
9:  c 
10:  c 

111  C 

12! 

13! 

1*1 

158 

16: 

I/: 

1  a : 
19: 
20: 
211 
22: 
23: 
24! 
251 
268 
278 
2s: 
29! 
30: 
31! 
32: 
33: 
3*: 

3b: 

36: 

37: 

3»: 
39:  1 
401 
*1: 

42: 

43:  2 
44: 

*51 
46!  3 

4/: 

48! 

49*. 

50! 

51 1 
52!  4 
53: 
54! 


PROGRAM  FOR  INITIAL  PRESSURE  DISTRIBUTION 
READ  ( 1  Ob* 14)  5#ALPRA»X«*XSTFP,iACTj*EFF#lDI«EN 
G»CiAMM*  (RATIO  OF  SPECIFIC  HEATS)  ' 

ALPHA  IS  di versing  NOZZLE  WALL  HALF  ANGLE 

XH  IS  LONGITUDINAL  RADIUS  OF  CURVATURE  pF  C^NVFRGIMU  SFCTION 
XSTfP  IS  STEP  SIZE  ALONG  X  COORDINATE  <PERCrNT  OF  L^CAL  -ALL  RAD) 
AACT  IS  THE  ACTUAL  NOZZLE  AREA  RATIO 

aeff  is  the  effectivf  nozzle  arfa  ratio 

IOIHEN-0  Fsr  T  ^0-DIHENSISNAL  NOZZLE  AND  1  F^R  AXISY^ETRJC  NOZZLE 
XK I G«  Z*IG,  AND  RWIG  AR£  THE  COORDINATES  OF  THE  PO^T  WHERE  THE 

NOZZLE  KALL  SLOPES  ARF  HATCHED  JUST  POWNSTRFAM  OF  MF  THROAT _ 

ALPHA* ALPHA *3. 14 159/ 130« 
nun*o 

PI2*3. 1*159/2. 

CBSALF*COSF (ALPHA) 

SINALF*SINF (ALPHA) 

.  I*aasjffl£lli*s  I N  ALEJ _ 

XwlG*XR* (PI ?♦ ALPHA ) 

RwIG*1»*XR*(1.-C8SAlF) 

xthoat*xr*pi? 

REND*SQRTF( AACT) 

IF  ( IDI",EN.tD»0)  RENO*aacT 

XEND«XV.IG-»(HENP.RW13)/SINALF _ 

C0NST*(  SSRfF  ( AACT )  -'SGR TF  { AEFF ) ) / ( COSALF*  ( Xr  M)Vxw  I G )  *  •  1 . 5 ) 

IF  (  lUlMEN.EU.O)  C9\ST*:AACT-AEFF)/(C0SALF*(XEND-XT1-,HAT)  ) 

ALPHA "ALPHA  *  180 ./ 3. 14159 
IF  (IDIHEN.EU.O)  WHITE  dOBf  13) 

IF  ( IUIhen.EU.I)  WRJTEI 1C8# 19) 

WRITE  (108.15)  3.ALPHA#XR,XENDfXSTEP. AACT, AEFF 
"  Cl  "•  ( H*l  *  ) /2  •  " 

C2*(G-10/P. 

C3*C1/(G-1.  ) 

C4*C1**C3 

C5*C4*C2*C3*2» 

FH.p.b _ _ 

X*O.Q 

R*0.0 

DS*0*0 

CONTINUE 

XU*X 

X*X*XSTtP*P 
IF  (X-Xw'lG)  2.2.4 
CONTINUE  ' 

IF  (  .NOT  •  (  X'j»LT«XTH0  AT  ,  AND.  X«GT ,  XTHO  AT  )  )  GO  TO  4 

X*XTH8A  r 

CONTINUE 

Z*XRM1.0-CBSF(X/yR> ) 

R«  1 , +XR*7  1  .-SINF(V/Wn - 

RINV*R 

C9SA*SINr(X/X«) 

GO  TO  5 

Z*(X-XWIG)«C9SALF+ZwIG 

R*(X-XWIG)»SINALF+RWIU 
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55: 

56: 

57! 

5h: 

5 

591 

63 1 

6 

611 

621 

63: 

6*: 

66: 

7 

66: 

8 

67: 

68! 

69: 

9 

70: 

7i: 

10 

72: 

11 

7i: 

7a: 

1? 

75: 

76: 

77: 

7a  1 

79: 

so: 

si : 

13 

82: 

63: 

C 

hi: 

1* 

85: 

15 

86: 

87: 

88 : 

89: 

16 

90: 

17 

911 

18 

921 

93: 

9ft 

19 

95: 

961 

971 

3S« CONST  *(X»>V»IG)**l1b  _  _ _  _  . . . . . 

IF  (IDI  "£"»•((*•  3)  L>S«C9NST«<X-X*IIG) 

RINV»M*05*C!JSALF 
A»{l«/«?IVV)**(l*Il>I  -IEN) 

IF  ( A-l -0)  7*6*6 
A-1.0 

FM1«I.0  ...  .  _ _ _ _ _ 

FM»1*0 

bM.O*C? 

08  T8  12 

IF  <X-PI2*XH)  11*11,9 
B“ ( 1 •♦C2*F 1  **2  1 

FM«FH^(A*n»»(C3»l>l»L»»FMl«  .)/<  CS»FM1»«2-C»»  _ 

IF  (AbSF|FM-FKH-0«00001)  12*11.11 
CONTINUE 

IF  (FK-1.20)  10*11*11 

FM«l*20 

FM1»FM 

GB  TU  8  _ _ 

CONTINUE 

PBP-U*»(G/(1*-G) ) 

num«num*i 

2«Z-XR 

WHITE  <108,161  FM,PAP,Z.X*D3,t<»8IKV.MUM 
WHITE  <106*171  X*H8P*MUH 
IF  (X-XEMD)  1*1*13' 

CONTINUE 

STOP 

FORMAT <6E10»0* 1 1 1 

F8HMAT(10X* 'GAMMA  ■ ' * F 7 . 3/* 1  OX* 

2 •  alpha  ■  ',F7 .?/.  ic nr; i m- - ■  r#T7'.'i/»iox.  'XE'^d — iT .rri'T/.  iox, 'XStep 

3  *'*F7.a/*10X*  *AACT  ■»»F7.1/»10X*'AEFF  •»*F7.3///«  8X*'M'*11X, 

4«P/P8'»  12X*  *Z'  *  13X*  •  X'  »12X.  'OS'.iaX,  'RW  *  U  X,  *RINV  *  7X,  'NUK'/l 
FORMAT  <IP7F14*5, 16) 

F0R-AT12E12.6,  53X,  13) 

FBRmAT ( 1  HI*  RX*  '8^E-DImEMSI9\AL  PEMFErT  QAS  EXPANSION  pHR  A  TWB» 

1  DIMENSIONAL  NOZZLE" GEOMETRY  a'TTW'VTTOX*  *  THF"A'5SUNPT  TON  THAT  THE  BO 

2Undary.laye«  displacement  t-mckness  varies  as  <x.xwtGit//i 
F8RMAT(1H1,RX* 'QnE-OIMEnSIP JAL  PERFECT  GAS  ; XPA-SIO1'  FOR  A%  AXISYM 
1  METRIC  NOZZLE  0EBMET&7  WITH  the*/,  iOX* ' ASSUMPT ION  THAT  THE  00 
2UNDAHY.LAYER  DISPLACEMENT  THICKNESS  VARIES  AS  <X-XWlG)»*3/g*//) 

END 
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